Merge pull request mom-ocean#51 from jiandewang/feature/update-to-GFD…

…L-20210120

Feature/update to gfdl 20210120

config_src/coupled_driver/MOM_surface_forcing_gfdl.F90

@@ -189,20 +189,13 @@ module MOM_surface_forcing_gfdl
                                                             !! ice-shelves, expressed as a coefficient
                                                             !! for divergence damping, as determined
                                                             !! outside of the ocean model [m3 s-1]
-  real, pointer, dimension(:,:)   :: ustk0           => NULL() !<
-  real, pointer, dimension(:,:)   :: vstk0           => NULL() !<
-  real, pointer, dimension(:)     :: stk_wavenumbers => NULL() !<
-  real, pointer, dimension(:,:,:) :: ustkb           => NULL() !<
-  real, pointer, dimension(:,:,:) :: vstkb           => NULL() !<
-
   integer :: xtype                    !< The type of the exchange - REGRID, REDIST or DIRECT
   type(coupler_2d_bc_type) :: fluxes  !< A structure that may contain an array of named fields
                                       !! used for passive tracer fluxes.
   integer :: wind_stagger = -999      !< A flag indicating the spatial discretization of wind stresses.
                                       !! This flag may be set by the flux-exchange code, based on what
                                       !! the sea-ice model is providing.  Otherwise, the value from
                                       !! the surface_forcing_CS is used.
-  integer :: num_stk_bands            !< Number of Stokes drift bands passed through the coupler
 end type ice_ocean_boundary_type
 
 integer :: id_clock_forcing !< A CPU time clock
@@ -684,7 +677,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
   real :: mass_eff      ! effective mass of sea ice for rigidity [R Z ~> kg m-2]
   real :: wt1, wt2      ! Relative weights of previous and current values of ustar, ND.
 
-  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0, istk
+  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer
   integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd
 
@@ -725,9 +718,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
        (associated(IOB%mass_berg) .and. (.not. associated(forces%mass_berg))) ) &
     call allocate_mech_forcing(G, forces, iceberg=.true.)
 
-  if ( associated(IOB%ustk0) ) &
-    call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=IOB%num_stk_bands)
-
   if (associated(IOB%ice_rigidity)) then
     rigidity_at_h(:,:) = 0.0
     call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed)
@@ -787,19 +777,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
       forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j)
     enddo ; enddo
   endif
-  forces%stk_wavenumbers(:) = IOB%stk_wavenumbers
-  do j=js,je; do i=is,ie
-     forces%ustk0(i,j) = IOB%ustk0(i-I0,j-J0) ! How to be careful here that the domains are right?
-     forces%vstk0(i,j) = IOB%vstk0(i-I0,j-J0)
-  enddo ; enddo
-  call pass_vector(forces%ustk0,forces%vstk0, G%domain )
-  do istk = 1,IOB%num_stk_bands
-    do j=js,je; do i=is,ie
-      forces%ustkb(i,j,istk) = IOB%ustkb(i-I0,j-J0,istk)
-      forces%vstkb(i,j,istk) = IOB%vstkb(i-I0,j-J0,istk)
-    enddo; enddo
-    call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
-  enddo
 
   ! Find the net mass source in the input forcing without other adjustments.
   if (CS%approx_net_mass_src .and. associated(forces%net_mass_src)) then

config_src/coupled_driver/ocean_model_MOM.F90

@@ -565,7 +565,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
     ! For now, the waves are only updated on the thermodynamics steps, because that is where
     ! the wave intensities are actually used to drive mixing.  At some point, the wave updates
     ! might also need to become a part of the ocean dynamics, according to B. Reichl.
-    call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces)
+    call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves)
   endif
 
   if ((OS%nstep==0) .and. (OS%nstep_thermo==0)) then ! This is the first call to update_ocean_model.

config_src/mct_driver/mom_surface_forcing_mct.F90

@@ -486,17 +486,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
     ! latent heat flux (W/m^2)
     fluxes%latent(i,j) = 0.0
-    ! contribution from frozen ppt
+    ! contribution from frozen ppt (notice minus sign since fprec is positive into the ocean)
     if (associated(IOB%fprec)) then
-      fluxes%latent(i,j)              = fluxes%latent(i,j) + &
+      fluxes%latent(i,j)              = fluxes%latent(i,j) - &
           IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
-      fluxes%latent_fprec_diag(i,j)   = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
+      fluxes%latent_fprec_diag(i,j)   = - G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
     endif
-    ! contribution from frozen runoff
+    ! contribution from frozen runoff (notice minus sign since rofi_flux is positive into the ocean)
     if (associated(fluxes%frunoff)) then
-      fluxes%latent(i,j)              = fluxes%latent(i,j) + &
+      fluxes%latent(i,j)              = fluxes%latent(i,j) - &
           IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
-      fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
+      fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
     endif
     ! contribution from evaporation
     if (associated(IOB%q_flux)) then

config_src/nuopc_driver/mom_cap.F90

@@ -36,6 +36,7 @@ module MOM_cap_mod
 use MOM_ocean_model_nuopc,    only: get_ocean_grid, get_eps_omesh
 use MOM_cap_time,             only: AlarmInit
 use MOM_cap_methods,          only: mom_import, mom_export, mom_set_geomtype, state_diagnose
+use MOM_cap_methods,          only: ChkErr
 #ifdef CESMCOUPLED
 use shr_file_mod,             only: shr_file_setLogUnit, shr_file_getLogUnit
 #endif
@@ -2049,18 +2050,6 @@ subroutine shr_file_getLogUnit(nunit)
 end subroutine shr_file_getLogUnit
 #endif
 
-logical function chkerr(rc, line, file)
-  integer, intent(in) :: rc
-  integer, intent(in) :: line
-  character(len=*), intent(in) :: file
-  integer :: lrc
-  chkerr = .false.
-  lrc = rc
-  if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then
-     chkerr = .true.
-  endif
-end function chkerr
-
 !>
 !! @page nuopc_cap NUOPC Cap
 !! @author Fei Liu ([email protected])

config_src/nuopc_driver/mom_cap_methods.F90

@@ -30,6 +30,7 @@ module MOM_cap_methods
 public :: mom_import
 public :: mom_export
 public :: state_diagnose
+public :: ChkErr
 
 private :: State_getImport
 private :: State_setExport
@@ -251,9 +252,9 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
 
 
   !----
-  ! Partitioned Stokes Drift Components 
+  ! Partitioned Stokes Drift Components
   !----
-  if ( associated(ice_ocean_boundary%ustkb) ) then    
+  if ( associated(ice_ocean_boundary%ustkb) ) then
     allocate(stkx1(isc:iec,jsc:jec))
     allocate(stky1(isc:iec,jsc:jec))
     allocate(stkx2(isc:iec,jsc:jec))
@@ -765,15 +766,18 @@ subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid
 
 end subroutine State_SetExport
 
+!> This subroutine writes the minimum, maximum and sum of each field
+!! contained within an ESMF state.
 subroutine state_diagnose(State, string, rc)
 
   ! ----------------------------------------------
   ! Diagnose status of State
   ! ----------------------------------------------
 
-  type(ESMF_State), intent(in)  :: state
-  character(len=*), intent(in)  :: string
-  integer         , intent(out) :: rc
+  type(ESMF_State), intent(in)    :: state  !< An ESMF State
+  character(len=*), intent(in)    :: string !< A string indicating whether the State is an
+                                            !! import or export State
+  integer         , intent(out)   :: rc     !< Return code
 
   ! local variables
   integer                         :: i,j,n
@@ -787,19 +791,19 @@ subroutine state_diagnose(State, string, rc)
   ! ----------------------------------------------
 
   call ESMF_StateGet(state, itemCount=fieldCount, rc=rc)
-  if (chkerr(rc,__LINE__,u_FILE_u)) return
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
   allocate(lfieldnamelist(fieldCount))
 
   call ESMF_StateGet(state, itemNameList=lfieldnamelist, rc=rc)
-  if (chkerr(rc,__LINE__,u_FILE_u)) return
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
   do n = 1, fieldCount
 
      call ESMF_StateGet(state, itemName=lfieldnamelist(n), field=lfield, rc=rc)
-     if (chkerr(rc,__LINE__,u_FILE_u)) return
+     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
      call field_getfldptr(lfield, fldptr1=dataPtr1d, fldptr2=dataPtr2d, rank=lrank, rc=rc)
-     if (chkerr(rc,__LINE__,u_FILE_u)) return
+     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
      if (lrank == 0) then
         ! no local data
@@ -829,23 +833,16 @@ subroutine state_diagnose(State, string, rc)
 
 end subroutine state_diagnose
 
-!===============================================================================
-
+!> Obtain a pointer to a rank 1 or rank 2 ESMF field
 subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
 
-  ! ----------------------------------------------
-  ! for a field, determine rank and return fldptr1 or fldptr2
-  ! abort is true by default and will abort if fldptr is not yet allocated in field
-  ! rank returns 0, 1, or 2.  0 means fldptr not allocated and abort=false
-  ! ----------------------------------------------
-
   ! input/output variables
-  type(ESMF_Field)  , intent(in)              :: field
-  real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr1(:)
-  real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr2(:,:)
-  integer           , intent(out)  , optional :: rank
-  logical           , intent(in)   , optional :: abort
-  integer           , intent(out)  , optional :: rc
+  type(ESMF_Field)  , intent(in)                        :: field        !< An ESMF field
+  real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr1(:)   !< A pointer to a rank 1 ESMF field
+  real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr2(:,:) !< A pointer to a rank 2 ESMF field
+  integer           , intent(out)            , optional :: rank         !< Field rank
+  logical           , intent(in)             , optional :: abort        !< Abort code
+  integer           , intent(out)            , optional :: rc           !< Return code
 
   ! local variables
   type(ESMF_GeomType_Flag)    :: geomtype
@@ -872,7 +869,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
   lrank = -99
 
   call ESMF_FieldGet(field, status=status, rc=rc)
-  if (chkerr(rc,__LINE__,u_FILE_u)) return
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
   if (status /= ESMF_FIELDSTATUS_COMPLETE) then
      lrank = 0
@@ -886,20 +883,20 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
   else
 
      call ESMF_FieldGet(field, geomtype=geomtype, rc=rc)
-     if (chkerr(rc,__LINE__,u_FILE_u)) return
+     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
      if (geomtype == ESMF_GEOMTYPE_GRID) then
         call ESMF_FieldGet(field, rank=lrank, rc=rc)
-        if (chkerr(rc,__LINE__,u_FILE_u)) return
+        if (ChkErr(rc,__LINE__,u_FILE_u)) return
      elseif (geomtype == ESMF_GEOMTYPE_MESH) then
         call ESMF_FieldGet(field, rank=lrank, rc=rc)
-        if (chkerr(rc,__LINE__,u_FILE_u)) return
+        if (ChkErr(rc,__LINE__,u_FILE_u)) return
         call ESMF_FieldGet(field, mesh=lmesh, rc=rc)
-        if (chkerr(rc,__LINE__,u_FILE_u)) return
+        if (ChkErr(rc,__LINE__,u_FILE_u)) return
         call ESMF_MeshGet(lmesh, numOwnedNodes=nnodes, numOwnedElements=nelements, rc=rc)
-        if (chkerr(rc,__LINE__,u_FILE_u)) return
+        if (ChkErr(rc,__LINE__,u_FILE_u)) return
         if (nnodes == 0 .and. nelements == 0) lrank = 0
-     else  
+     else
         call ESMF_LogWrite(trim(subname)//": ERROR geomtype not supported ", &
              ESMF_LOGMSG_INFO)
         rc = ESMF_FAILURE
@@ -917,7 +914,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
            return
         endif
         call ESMF_FieldGet(field, farrayPtr=fldptr1, rc=rc)
-        if (chkerr(rc,__LINE__,u_FILE_u)) return
+        if (ChkErr(rc,__LINE__,u_FILE_u)) return
      elseif (lrank == 2) then
         if (.not.present(fldptr2)) then
            call ESMF_LogWrite(trim(subname)//": ERROR missing rank=2 array ", &
@@ -926,7 +923,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
            return
         endif
         call ESMF_FieldGet(field, farrayPtr=fldptr2, rc=rc)
-        if (chkerr(rc,__LINE__,u_FILE_u)) return
+        if (ChkErr(rc,__LINE__,u_FILE_u)) return
      else
         call ESMF_LogWrite(trim(subname)//": ERROR in rank ", &
              ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
@@ -942,16 +939,17 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
 
 end subroutine field_getfldptr
 
-logical function chkerr(rc, line, file)
-  integer, intent(in) :: rc
-  integer, intent(in) :: line
-  character(len=*), intent(in) :: file
+!> Returns true if ESMF_LogFoundError() determines that rc is an error code. Otherwise false.
+logical function ChkErr(rc, line, file)
+  integer, intent(in) :: rc            !< return code to check
+  integer, intent(in) :: line          !< Integer source line number
+  character(len=*), intent(in) :: file !< User-provided source file name
   integer :: lrc
-  chkerr = .false.
+  ChkErr = .false.
   lrc = rc
   if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then
-     chkerr = .true.
+     ChkErr = .true.
   endif
-end function chkerr
+end function ChkErr
 
 end module MOM_cap_methods

config_src/nuopc_driver/mom_cap_time.F90

@@ -16,6 +16,7 @@ module MOM_cap_time
 use ESMF                  , only : ESMF_RC_ARG_BAD
 use ESMF                  , only : operator(<), operator(/=), operator(+), operator(-), operator(*) , operator(>=)
 use ESMF                  , only : operator(<=), operator(>), operator(==)
+use MOM_cap_methods       , only : ChkErr
 
 implicit none; private
 
@@ -336,16 +337,4 @@ subroutine date2ymd (date, year, month, day)
 
 end subroutine date2ymd
 
-logical function chkerr(rc, line, file)
-  integer, intent(in) :: rc
-  integer, intent(in) :: line
-  character(len=*), intent(in) :: file
-  integer :: lrc
-  chkerr = .false.
-  lrc = rc
-  if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then
-     chkerr = .true.
-  endif
-end function chkerr
-
 end module MOM_cap_time

config_src/nuopc_driver/mom_surface_forcing_nuopc.F90

@@ -183,11 +183,15 @@ module MOM_surface_forcing_nuopc
                                                               !! ice-shelves, expressed as a coefficient
                                                               !! for divergence damping, as determined
                                                               !! outside of the ocean model in [m3/s]
-  real, pointer, dimension(:,:)   :: ustk0           => NULL() !<
-  real, pointer, dimension(:,:)   :: vstk0           => NULL() !<
-  real, pointer, dimension(:)     :: stk_wavenumbers => NULL() !<
-  real, pointer, dimension(:,:,:) :: ustkb           => NULL() !<
-  real, pointer, dimension(:,:,:) :: vstkb           => NULL() !<
+  real, pointer, dimension(:,:)   :: ustk0           => NULL() !< Surface Stokes drift, zonal [m/s]
+  real, pointer, dimension(:,:)   :: vstk0           => NULL() !< Surface Stokes drift, meridional [m/s]
+  real, pointer, dimension(:)     :: stk_wavenumbers => NULL() !< The central wave number of Stokes bands [rad/m]
+  real, pointer, dimension(:,:,:) :: ustkb           => NULL() !< Stokes Drift spectrum, zonal [m/s]
+                                                               !! Horizontal  - u points
+                                                               !! 3rd dimension - wavenumber
+  real, pointer, dimension(:,:,:) :: vstkb           => NULL() !< Stokes Drift spectrum, meridional [m/s]
+                                                               !! Horizontal  - v points
+                                                               !! 3rd dimension - wavenumber
   integer :: num_stk_bands            !< Number of Stokes drift bands passed through the coupler
   integer :: xtype                                            !< The type of the exchange - REGRID, REDIST or DIRECT
   type(coupler_2d_bc_type)      :: fluxes                     !< A structure that may contain an array of
@@ -493,15 +497,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
          fluxes%seaice_melt(i,j) = kg_m2_s_conversion * G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0)
 
     fluxes%latent(i,j) = 0.0
+    ! notice minus sign since fprec is positive into the ocean
     if (associated(IOB%fprec)) then
-       fluxes%latent(i,j)            = fluxes%latent(i,j) + &
+       fluxes%latent(i,j)            = fluxes%latent(i,j) - &
            IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
-       fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
+       fluxes%latent_fprec_diag(i,j) = - G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
     endif
+    ! notice minus sign since frunoff is positive into the ocean
     if (associated(IOB%frunoff)) then
-       fluxes%latent(i,j)              = fluxes%latent(i,j) + &
+       fluxes%latent(i,j)              = fluxes%latent(i,j) - &
            IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
-       fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * &
+       fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * &
            IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
     endif
     if (associated(IOB%q_flux)) then
@@ -858,7 +864,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
       enddo; enddo
       call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
     enddo
-  endif 
+  endif
 
   ! sea ice related dynamic fields
   if (associated(IOB%ice_rigidity)) then

config_src/nuopc_driver/time_utils.F90

@@ -14,6 +14,7 @@ module time_utils_mod
 use ESMF,               only: ESMF_Time, ESMF_TimeGet, ESMF_LogFoundError
 use ESMF,               only: ESMF_LOGERR_PASSTHRU,ESMF_TimeInterval
 use ESMF,               only: ESMF_TimeIntervalGet, ESMF_TimeSet, ESMF_SUCCESS
+use MOM_cap_methods,    only: ChkErr
 
 implicit none; private
 
@@ -160,16 +161,4 @@ function string_to_date(string, rc)
 
 end function string_to_date
 
-logical function chkerr(rc, line, file)
-  integer, intent(in) :: rc
-  integer, intent(in) :: line
-  character(len=*), intent(in) :: file
-  integer :: lrc
-  chkerr = .false.
-  lrc = rc
-  if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then
-     chkerr = .true.
-  endif
-end function chkerr
-
 end module time_utils_mod