NOAA-GFDL / SHiELD_physics

Driver for the SHiELD model
Other
3 stars 11 forks source link

how to produce diagnostic winds at time step 0? #32

Open StevePny opened 11 months ago

StevePny commented 11 months ago

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 11 months ago

Hi, Steve. Thank you for opening this issue. I will take a closer look upon returning from travel but I believe the upshot is that u_srf/v_srf need to have the lat-lon winds to have been previously computed; this is done at the end of fv_dynamics() but not necessarily in the model initialization.

These are indeed the lowest-layer winds which are not the same as the 10-m winds; the latter is computed from the physics by Monin-Obhukov theory.

Lucas

On Mon, Nov 13, 2023 at 9:10 PM Steve Penny @.***> wrote:

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.

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/SHiELD_physics/issues/32, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVGQTBAF5QLMRLOCLTDYELHHTAVCNFSM6AAAAAA7KCHD32VHI2DSMVQWIX3LMV43ASLTON2WKOZRHE4TCOBUGM4DOMQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

StevePny commented 11 months ago

Thanks Lucas (@lharris4), we are currently using the 10m winds computed here, and they appear to be accurate (we have tested against scatterometer observations): https://github.com/NOAA-GFDL/SHiELD_physics/blob/a5a3a40c6fa656097e3de47a89a44e619effa399/gsmphys/sfc_diag.f#L37-L38

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: https://github.com/NOAA-GFDL/SHiELD_physics/blob/2882fdeb429abc2349a8e881803ac67b154532c3/simple_coupler/coupler_main.F90#L216

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 11 months ago

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 11 months ago

Hi, Steve. Thank you for sending this information. I am juggling a few things during this holiday week but hope to get back to this shortly.

Lucas

On Fri, Nov 17, 2023 at 8:16 PM Steve Penny @.***> wrote:

FYI @lharris4 https://github.com/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: [image: Screen Shot 2023-11-17 at 3 52 11 PM] https://user-images.githubusercontent.com/11303168/283951346-9daa1628-c75f-4f9e-ae0b-af753abadba0.png [image: Screen Shot 2023-11-17 at 3 52 43 PM] https://user-images.githubusercontent.com/11303168/283951359-d7a14ae4-0513-454c-a946-95e019550c90.png

— Reply to this email directly, view it on GitHub https://github.com/NOAA-GFDL/SHiELD_physics/issues/32#issuecomment-1817306244, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMUQRVFMSWYHMUNBYUYHMXLYFAD73AVCNFSM6AAAAAA7KCHD32VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMJXGMYDMMRUGQ . You are receiving this because you were mentioned.Message ID: @.***>

lharris4 commented 11 months ago

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 11 months ago

3) The missing u_srf appears to be a data synchronization issue. Could you show me the code that produces this output file?

StevePny commented 10 months ago

@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 8 months ago

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