NOAA-GFDL/SHiELD_physics

how to produce diagnostic winds at time step 0?

Opened this issue · 8 comments

StevePny commented

Is your question related to a problem? Please describe.
From our own external driver routine, we are accessing the 10m winds through IPD_Data(nb)%intdiag%u10m and IPD_Data(nb)%intdiag%v10m. This works fine for us after the first model time step. However we have noticed that on time step 0 (upon initialization of the model) these fields are empty. It appears that they are diagnostic fields that only get populated after the first model time step completes. We would like realistic u10m/v10m fields to be populated at time 0, e.g. based on the initial condition file or restart file.

Is there a straightforward way to trigger the diagnostic computation of u10m/v10m winds when the model initializes so that they are not empty at time step 0?

Describe what you have tried
We've tried accessing the Atm(mygrid)%u_srf and Atm(mygrid)%v_srf fields from the fv3 dynamical core at timestep 0, but they appear to be represented slightly differently on the grid. It's not clear to me whether these are on a C-grid or D-grid stagger while it appears the diagnostic 10m winds above are at cell centers. Ideally, we'd like them to be on the same grid at the same associated height. The bigger problem for us in this case was identifying the appropriate means of collecting the parallel-distributed data, which appears to be different from what was used above for the IPD_Data(nb)%intdiag data structure.

We attempted something like:

      if (varchar=='u') then
          do j=jsc,jec
              do i=isc,iec
                  dataPtr_r8(i,j) = real(Atm(mygrid)%u_srf(i,j),kind=8)
              enddo
          enddo
      elseif (varchar=='v') then
          do j=jsc,jec
              do i=isc,iec
                  dataPtr_r8(i,j) = real(Atm(mygrid)%v_srf(i,j),kind=8)
              enddo
          enddo
      endif

but the resulting fields were not reconstructed correctly. Further, we presume these fields are defined at the model's bottom level and not at the desired 10m height.

lharris4 commented
StevePny commented

Thanks Lucas (@lharris4), we are currently using the 10m winds computed here, and they appear to be accurate (we have tested against scatterometer observations):

SHiELD_physics/gsmphys/sfc_diag.f

Lines 37 to 38 in a5a3a40

u10m(i) = f10m(i) * u1(i)
v10m(i) = f10m(i) * v1(i)

For example, here is tile 1 after 1 model time step:
Screen Shot 2023-11-16 at 5 30 37 PM

and what the wave model is seeing at the same time step:
Screen Shot 2023-11-16 at 5 31 53 PM

Our driver routine runs this initialization sequence, which was copied from the SHiELD "program coupler_main":

    call fms_init(mpi_comm_fv3)
    call mpp_init()
    initClock = mpp_clock_id( 'Initialization' )
    call mpp_clock_begin (initClock)

    call constants_init
    call sat_vapor_pres_init

    call coupler_init(model=model,clock=dclock,rc=rc)

where 'model' and 'dclock' are ESMF objects, and our 'coupler_init' is as close as possible to:

SHiELD_physics/simple_coupler/coupler_main.F90

Line 216 in 2882fde

subroutine coupler_init

At initialization time, it appears that all data stored in "IPD_Data(nb)%Sfcprop" have been populated with values, while all data stored in IPD_Data(nb)%intdiag are empty, hence our challenge with 10m winds at time step 0.

StevePny commented

FYI @lharris4 this first plot is what we see when accessing Atm(mygrid)%u_srf after initialization, using the loops described in the original post above. The second plot shows u10m diagnostic (intdiag) at time 1 (it is empty at initialization time 0). I'm not sure why we're losing half the field, but if we can recover it then we can replicate the u10m/v10m diagnostic Monin-Obhukov computation to get what we need:

Screen Shot 2023-11-17 at 3 52 11 PM Screen Shot 2023-11-17 at 3 52 43 PM
lharris4 commented
lharris4 commented

Hi @StevePny. I think there are a couple of convolved issues here.

  1. u_srf/v_srf are indeed lowest-layer and not 10-m winds, although 10-m winds can be computed using M-O theory from the lowest-layer winds. It is a bit tricky at initialization since you would need to compute f10m yourself (it may not yet be initialized in the physics), but doable. The diagnostics from the physics will indeed not be available at initialization unless you add extra code to the SHiELD physics initialization to do this.

  2. Atm(mygrid)%u_srf should be available after calling fv_restart(), as it will either be copied from Atm%u or read from the restart files. This will be the zonal (earth-relative) wind defined as a cell-centered value, and not staggered or in a grid-relative coordinate.

(next comment to follow)

lharris4 commented
  1. The missing u_srf appears to be a data synchronization issue. Could you show me the code that produces this output file?
StevePny commented

@lharris4
Sure - we created a simple output function here. We are just copying the functionality of "get_bottom_wind ( u_bot, v_bot )", but with the option to choose 'u' or 'v' independently:

    subroutine get_srfwinds(dataPtr_r8, varchar)
!     use GFS_typedefs,       only: kind_phys
      use atmos_model_mod, only: Atm_block
      use atmosphere_mod, only: Atm, mygrid, get_bottom_wind

      real(ESMF_KIND_R8), dimension(:,:), pointer, intent(inout) :: dataPtr_r8
      character(1), intent(in) :: varchar  ! either 'u' or 'v'
      integer :: isc,iec,jsc,jec
      integer :: i,j
!     real(ESMF_KIND_R4), dimension(:,:), allocatable :: u_bot, v_bot

      logical :: dodebug = .false.

      ! -------------------------------------------------------------------------
      ! Check fv3 dynamical core arrays to make sure they are available
      ! -------------------------------------------------------------------------
      if (.not. allocated(Atm(mygrid)%u_srf) ) then
          print *, "fv3_shield_cap::setExport::get_srfwinds: ERROR - Atm(mygrid)%u_srf is not allocated."
          print *, "EXITING..."
          stop(1)
      elseif (.not. allocated(Atm(mygrid)%v_srf) ) then
          print *, "fv3_shield_cap::setExport::get_srfwinds: ERROR - Atm(mygrid)%v_srf is not allocated."
          print *, "EXITING..."
          stop(1)
      endif

      ! -------------------------------------------------------------------------
      ! Get array indices
      ! -------------------------------------------------------------------------
      isc = Atm_block%isc
      iec = Atm_block%iec
      jsc = Atm_block%jsc
      jec = Atm_block%jec
  !   nk  = Atm_block%npz

      ! -------------------------------------------------------------------------
      ! Assign the atmos fields to the ESMF data pointers
      ! -----------------------
      !STEVE:NOTE:: https://github.com/NOAA-EMC/fv3atm/blob/deeac5f0acb875f8a282200ebdc69d6157232163/atmos_model.F90#L2829
      ! -----------------------
      ! Set the u/v winds fields 
      ! using Atm(mygrid)%u_srf and Atm(mygrid)%v_srf
      !
      ! Note: the ATM data structure comes from the fv3 dynamical core:
      ! https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/blob/efcb02f06dec753518e7fb51c90e300d36a7455f/driver/SHiELD/atmosphere.F90#L1112C13-L1112C45
      ! -------------------------------------------------------------------------
!     if (.not.allocated(u_bot)) allocate(u_bot(isc:iec,jsc:jec))
!     if (.not.allocated(v_bot)) allocate(v_bot(isc:iec,jsc:jec))

!     call get_bottom_wind ( u_bot, v_bot )
      if (varchar=='u') then
          do j=jsc,jec
              do i=isc,iec
                  dataPtr_r8(i,j) = real(Atm(mygrid)%u_srf(i,j),kind=8)
              enddo
          enddo
      elseif (varchar=='v') then
          do j=jsc,jec
              do i=isc,iec
                  dataPtr_r8(i,j) = real(Atm(mygrid)%v_srf(i,j),kind=8)
              enddo
          enddo
      endif

!     if (allocated(u_bot)) deallocate(u_bot)
!     if (allocated(v_bot)) deallocate(v_bot)

    end subroutine get_srfwinds

Which is written to file with:

        ! -----------------------
        ! Write all fields in the export state
        ! -----------------------
        call NUOPC_Write(exportState, fileNamePrefix="ATM_DataInitialize_exportState_tile*_", overwrite=.true., rc=rc)
        if (CheckError(rc,__LINE__,__FILE__)) return
lharris4 commented

Hi @StevePny . Sorry for the delay. I don't see anything in this code that might cause a problem. The only things that come to my mind are (a) that a sync didn't complete by the time of the write, or (b) for some reason the srf_init flag could sometimes act in a funny way. Have you been able to try this functionality again recently?

Thanks
Lucas