Dev master candidate from NCAR - 2019 06 11

config_src/mct_driver/MOM_surface_forcing.F90

@@ -215,7 +215,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
   type(forcing),           intent(inout) :: fluxes !< A structure containing pointers to
                                                    !! all possible mass, heat or salt flux forcing fields.
                                                    !!  Unused fields have NULL ptrs.
-
   type(time_type),         intent(in)    :: Time   !< The time of the fluxes, used for interpolating the
                                                    !! salinity to the right time, when it is being restored.
   type(ocean_grid_type),   intent(inout) :: G      !< The ocean's grid structure
@@ -244,7 +243,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
 
   integer :: i, j, k, 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
 
   logical :: restore_salinity ! local copy of the argument restore_salt, if it
                               ! is present, or false (no restoring) otherwise.
@@ -395,9 +393,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
     enddo; enddo
   endif
 
-  !i0 = is - isc_bnd ; j0 = js - jsc_bnd ???
-  i0 = 0; j0 = 0 ! TODO: is this right?
-
+  ! obtain fluxes from IOB
+  i0 = 0; j0 = 0
   do j=js,je ; do i=is,ie
     ! liquid precipitation (rain)
     if (associated(fluxes%lprec)) &
@@ -456,8 +453,27 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, US, CS, &
       fluxes%seaice_melt(i,j) = G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0)
 
     ! latent heat flux (W/m^2)
-    if (associated(fluxes%latent)) &
-      fluxes%latent(i,j) = G%mask2dT(i,j) * IOB%latent_flux(i-i0,j-j0)
+    ! old method, latent = IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+    !if (associated(fluxes%latent)) &
+    !  fluxes%latent(i,j) = G%mask2dT(i,j) * IOB%latent_flux(i-i0,j-j0)
+    ! new method
+    fluxes%latent(i,j) = 0.0
+    ! contribution from frozen ppt
+    if (associated(fluxes%fprec)) then
+      fluxes%latent(i,j)              = fluxes%latent(i,j) + IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion
+      fluxes%latent_fprec_diag(i,j)   = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion
+    endif
+    ! contribution from frozen runoff
+    if (associated(fluxes%frunoff)) then
+      fluxes%latent(i,j)              = fluxes%latent(i,j) + IOB%rofi_flux(i-i0,j-j0)*CS%latent_heat_fusion
+      fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*CS%latent_heat_fusion
+    endif
+    ! contribution from evaporation
+    if (associated(IOB%q_flux)) then
+      fluxes%latent(i,j)             = fluxes%latent(i,j) + IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+      fluxes%latent_evap_diag(i,j)  = G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+    endif
+    fluxes%latent(i,j) = G%mask2dT(i,j) * fluxes%latent(i,j)
 
     if (associated(IOB%sw_flux_vis_dir)) &
       fluxes%sw_vis_dir(i,j) = G%mask2dT(i,j) * IOB%sw_flux_vis_dir(i-i0,j-j0)
@@ -580,8 +596,9 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
 
  !isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2)
  !jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4)
+ !if (is_root_pe()) write(*,*)'isc_bnd, jsc_bnd, iec_bnd, jec_bnd',isc_bnd, jsc_bnd, iec_bnd, jec_bnd
  !i0 = is - isc_bnd ; j0 = js - jsc_bnd
-  i0 = 0; j0 = 0 ! TODO: is this right?
+ i0 = 0; j0 = 0 ! TODO: is this right?
 
   Irho0 = 1.0/CS%Rho0
 

config_src/mct_driver/ocn_cap_methods.F90

@@ -50,11 +50,11 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
       ! rotate taux and tauy from true zonal/meridional to local coordinates
       ! taux
       ice_ocean_boundary%u_flux(i,j) = GRID%cos_rot(i,j) * x2o(ind%x2o_Foxx_taux,k) &
-                                      + GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k)
+                                      - GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k)
 
       ! tauy
       ice_ocean_boundary%v_flux(i,j) = GRID%cos_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k) &
-                                      - GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_taux,k)
+                                      + GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_taux,k)
 
       ! liquid precipitation (rain)
       ice_ocean_boundary%lprec(i,j) = x2o(ind%x2o_Faxa_rain,k)
@@ -188,9 +188,9 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
       o2x(ind%o2x_So_t, n) = ocn_public%t_surf(ig,jg) * grid%mask2dT(i,j)
       o2x(ind%o2x_So_s, n) = ocn_public%s_surf(ig,jg) * grid%mask2dT(i,j)
       ! rotate ocn current from local tripolar grid to true zonal/meridional (inverse transformation)
-      o2x(ind%o2x_So_u, n) = (grid%cos_rot(i,j) * ocn_public%u_surf(ig,jg) - &
+      o2x(ind%o2x_So_u, n) = (grid%cos_rot(i,j) * ocn_public%u_surf(ig,jg) + &
                               grid%sin_rot(i,j) * ocn_public%v_surf(ig,jg)) * grid%mask2dT(i,j)
-      o2x(ind%o2x_So_v, n) = (grid%cos_rot(i,j) * ocn_public%v_surf(ig,jg) + &
+      o2x(ind%o2x_So_v, n) = (grid%cos_rot(i,j) * ocn_public%v_surf(ig,jg) - &
                               grid%sin_rot(i,j) * ocn_public%u_surf(ig,jg)) * grid%mask2dT(i,j)
 
       ! boundary layer depth (m)
@@ -270,8 +270,8 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
   n = 0
   do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
     n = n+1
-    o2x(ind%o2x_So_dhdx, n) = grid%cos_rot(i,j) * sshx(i,j) - grid%sin_rot(i,j) * sshy(i,j)
-    o2x(ind%o2x_So_dhdy, n) = grid%cos_rot(i,j) * sshy(i,j) + grid%sin_rot(i,j) * sshx(i,j)
+    o2x(ind%o2x_So_dhdx, n) = grid%cos_rot(i,j) * sshx(i,j) + grid%sin_rot(i,j) * sshy(i,j)
+    o2x(ind%o2x_So_dhdy, n) = grid%cos_rot(i,j) * sshy(i,j) - grid%sin_rot(i,j) * sshx(i,j)
   enddo; enddo
 
 end subroutine ocn_export

config_src/nuopc_driver/MOM_surface_forcing.F90

@@ -163,7 +163,7 @@ module MOM_surface_forcing
   real, pointer, dimension(:,:) :: q_flux            =>NULL() !< specific humidity flux [kg/m2/s]
   real, pointer, dimension(:,:) :: salt_flux         =>NULL() !< salt flux [kg/m2/s]
   real, pointer, dimension(:,:) :: seaice_melt_heat  =>NULL() !< sea ice and snow melt heat flux [W/m2]
-  real, pointer, dimension(:,:) :: seaice_melt_water =>NULL() !< water flux due to sea ice and snow melting [kg/m2/s]
+  real, pointer, dimension(:,:) :: seaice_melt       =>NULL() !< water flux due to sea ice and snow melting [kg/m2/s]
   real, pointer, dimension(:,:) :: lw_flux           =>NULL() !< long wave radiation [W/m2]
   real, pointer, dimension(:,:) :: sw_flux_vis_dir   =>NULL() !< direct visible sw radiation [W/m2]
   real, pointer, dimension(:,:) :: sw_flux_vis_dif   =>NULL() !< diffuse visible sw radiation [W/m2]
@@ -457,18 +457,18 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, &
          fluxes%heat_content_frunoff(i,j) = IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j)
 
     if (associated(IOB%lw_flux)) &
-      fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j)
+         fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j)
 
     if (associated(IOB%t_flux)) &
-      fluxes%sens(i,j) = IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j)
+         fluxes%sens(i,j) = IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j)
 
-    ! ! sea ice and snow melt heat flux [W/m2]
-    ! if (associated(fluxes%seaice_melt_heat)) &
-    !   fluxes%seaice_melt_heat(i,j) = G%mask2dT(i,j) * IOB%seaice_melt_heat(i-i0,j-j0)
+    ! sea ice and snow melt heat flux [W/m2]
+    if (associated(IOB%seaice_melt_heat)) &
+         fluxes%seaice_melt_heat(i,j) = G%mask2dT(i,j) * IOB%seaice_melt_heat(i-i0,j-j0)
 
-    ! ! water flux due to sea ice and snow melt [kg/m2/s]
-    ! if (associated(fluxes%seaice_melt)) &
-    !   fluxes%seaice_melt(i,j) = G%mask2dT(i,j) * IOB%seaice_melt_water(i-i0,j-j0)
+    ! water flux due to sea ice and snow melt [kg/m2/s]
+    if (associated(IOB%seaice_melt)) &
+         fluxes%seaice_melt(i,j) = G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0)
 
     fluxes%latent(i,j) = 0.0
     if (associated(IOB%fprec)) then
@@ -540,10 +540,9 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, US, CS, &
     sign_for_net_FW_bug = 1.
     if (CS%use_net_FW_adjustment_sign_bug) sign_for_net_FW_bug = -1.
     do j=js,je ; do i=is,ie
-      net_FW(i,j) = (((fluxes%lprec(i,j)   + fluxes%fprec(i,j)) + &
+      net_FW(i,j) = (((fluxes%lprec(i,j)   + fluxes%fprec(i,j) + fluxes%seaice_melt(i,j)) + &
                       (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
                       (fluxes%evap(i,j)    + fluxes%vprec(i,j)) ) * G%areaT(i,j)
-      ! net_FW(i,j) = netFW(i,j) + fluxes%seaice_melt(i,j) * G%areaT(i,j)
 
       !   The following contribution appears to be calculating the volume flux of sea-ice
       ! melt. This calculation is clearly WRONG if either sea-ice has variable
@@ -1087,7 +1086,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt,
     call get_param(param_file, mdl, "USE_NET_FW_ADJUSTMENT_SIGN_BUG", &
                  CS%use_net_FW_adjustment_sign_bug, &
                    "If true, use the wrong sign for the adjustment to\n"//&
-                   "the net fresh-water.", default=.true.)
+                   "the net fresh-water.", default=.false.)
   call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_BY_SCALING", &
                  CS%adjust_net_fresh_water_by_scaling, &
                  "If true, adjustments to net fresh water to achieve zero net are\n"//&
@@ -1368,8 +1367,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   write(outunit,100) 'iobt%t_flux         '   , mpp_chksum( iobt%t_flux         )
   write(outunit,100) 'iobt%q_flux         '   , mpp_chksum( iobt%q_flux         )
   write(outunit,100) 'iobt%salt_flux      '   , mpp_chksum( iobt%salt_flux      )
- !write(outunit,100) 'iobt%seaice_melt_heat'  , mpp_chksum( iobt%seaice_melt_heat)
- !write(outunit,100) 'iobt%seaice_melt_water' , mpp_chksum( iobt%seaice_melt_water)
+  write(outunit,100) 'iobt%seaice_melt_heat'  , mpp_chksum( iobt%seaice_melt_heat)
+  write(outunit,100) 'iobt%seaice_melt    '   , mpp_chksum( iobt%seaice_melt    )
   write(outunit,100) 'iobt%lw_flux        '   , mpp_chksum( iobt%lw_flux        )
   write(outunit,100) 'iobt%sw_flux_vis_dir'   , mpp_chksum( iobt%sw_flux_vis_dir)
   write(outunit,100) 'iobt%sw_flux_vis_dif'   , mpp_chksum( iobt%sw_flux_vis_dif)

config_src/nuopc_driver/mom_cap.F90

@@ -204,6 +204,8 @@
 !! --------------------------|------------|-----------------|---------------------------------------|-------------------
 !! inst_pres_height_surface  | Pa         | p               | pressure of overlying sea ice and atmosphere
 !! mass_of_overlying_sea_ice | kg         | mi              | mass of overlying sea ice          | |
+!! seaice_melt_heat          | W m-2      | seaice_melt_heat| sea ice and snow melt heat flux    | |
+!! seaice_melt               | kg m-2 s-1 | seaice_melt     | water flux due to sea ice and snow melting    | |
 !! mean_calving_heat_flx     | W m-2      | calving_hflx    | heat flux, relative to 0C, of frozen land water into ocean
 !! mean_calving_rate         | kg m-2 s-1 | calving         | mass flux of frozen runoff         | |
 !! mean_evap_rate            | kg m-2 s-1 | q_flux          | specific humidity flux             |
@@ -377,6 +379,7 @@ module mom_cap_mod
 use ESMF,  only: ESMF_COORDSYS_SPH_DEG, ESMF_GridCreate, ESMF_INDEX_DELOCAL
 use ESMF,  only: ESMF_MESHLOC_ELEMENT, ESMF_RC_VAL_OUTOFRANGE, ESMF_StateGet
 use ESMF,  only: ESMF_TimePrint, ESMF_AlarmSet, ESMF_FieldGet
+use ESMF,  only: operator(==), operator(/=), operator(+), operator(-)
 
 ! TODO ESMF_GridCompGetInternalState does not have an explicit Fortran interface.
 !! Model does not compile with "use ESMF,  only: ESMF_GridCompGetInternalState"
@@ -753,6 +756,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   integer                                :: userRc
   character(len=512)                     :: restartfile          ! Path/Name of restart file
   character(len=*), parameter            :: subname='(mom_cap:InitializeAdvertise)'
+  character(len=32)                      :: calendar
 !--------------------------------
 
   rc = ESMF_SUCCESS
@@ -804,7 +808,35 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fms_init(mpi_comm_mom)
   call constants_init
   call field_manager_init
-  call set_calendar_type (JULIAN)
+
+  ! determine the calendar
+  if (cesm_coupled) then
+     call NUOPC_CompAttributeGet(gcomp, name="calendar", value=cvalue, &
+          isPresent=isPresent, isSet=isSet, rc=rc)
+     if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+          line=__LINE__, &
+          file=__FILE__)) &
+          return  ! bail out
+     if (isPresent .and. isSet) then
+        read(cvalue,*) calendar
+        select case (trim(calendar))
+           case ("NO_LEAP")
+              call set_calendar_type (NOLEAP)
+           case ("GREGORIAN")
+              call set_calendar_type (GREGORIAN)
+           case default
+              call ESMF_LogSetError(ESMF_RC_ARG_BAD, &
+                 msg=subname//": Calendar not supported in MOM6: "//trim(calendar), &
+                 line=__LINE__, file=__FILE__, rcToReturn=rc)
+           end select
+     else
+        call set_calendar_type (NOLEAP)
+     endif
+
+  else
+     call set_calendar_type (JULIAN)
+  endif
+
   call diag_manager_init
 
   ! this ocean connector will be driven at set interval
@@ -960,6 +992,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
              Ice_ocean_boundary% sw_flux_nir_dif (isc:iec,jsc:jec), &
              Ice_ocean_boundary% lprec (isc:iec,jsc:jec),           &
              Ice_ocean_boundary% fprec (isc:iec,jsc:jec),           &
+             Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),&
+             Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec),     &
              Ice_ocean_boundary% mi (isc:iec,jsc:jec),              &
              Ice_ocean_boundary% p (isc:iec,jsc:jec),               &
              Ice_ocean_boundary% runoff (isc:iec,jsc:jec),          &
@@ -981,6 +1015,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   Ice_ocean_boundary%sw_flux_nir_dif = 0.0
   Ice_ocean_boundary%lprec           = 0.0
   Ice_ocean_boundary%fprec           = 0.0
+  Ice_ocean_boundary%seaice_melt     = 0.0
+  Ice_ocean_boundary%seaice_melt_heat= 0.0
   Ice_ocean_boundary%mi              = 0.0
   Ice_ocean_boundary%p               = 0.0
   Ice_ocean_boundary%runoff          = 0.0
@@ -1030,8 +1066,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "inst_pres_height_surface"   , "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl"                  , "will provide") !-> liquid runoff
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi"                  , "will provide") !-> ice runoff
-  !call fld_list_add(fldsToOcn_num, fldsToOcn, "seaice_melt_water"          , "will provide")
-  !call fld_list_add(fldsToOcn_num, fldsToOcn, "seaice_melt_heat"           , "will provide")
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide")
+  call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn"        , "will provide")
 
  !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_rate"           , "will provide")
  !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_rate"          , "will provide")
@@ -1728,6 +1764,7 @@ subroutine ModelAdvance(gcomp, rc)
   ! local variables
   integer                                :: userRc
   logical                                :: existflag, isPresent, isSet
+  logical                                :: do_advance = .true.
   type(ESMF_Clock)                       :: clock!< ESMF Clock class definition
   type(ESMF_Alarm)                       :: alarm
   type(ESMF_State)                       :: importState, exportState
@@ -1810,8 +1847,47 @@ subroutine ModelAdvance(gcomp, rc)
     file=__FILE__)) &
     return  ! bail out
 
-  Time = esmf2fms_time(currTime)
   Time_step_coupled = esmf2fms_time(timeStep)
+  Time = esmf2fms_time(currTime)
+
+  !---------------
+  ! Apply ocean lag for startup runs:
+  !---------------
+
+  if (cesm_coupled) then
+    if (trim(runtype) == "initial") then
+
+      ! Do not call MOM6 timestepping routine if the first cpl tstep of a startup run
+      if (currTime == startTime) then
+        call ESMF_LogWrite("MOM6 - Skipping the first coupling timestep", ESMF_LOGMSG_INFO, rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+          line=__LINE__, &
+          file=__FILE__)) &
+          return  ! bail out
+        do_advance = .false.
+      else
+        do_advance = .true.
+      endif
+
+      ! If the second cpl tstep of a startup run, step back a cpl tstep and advance for two cpl tsteps
+      if (currTime == startTime + timeStep) then
+        call ESMF_LogWrite("MOM6 - Stepping back one coupling timestep", ESMF_LOGMSG_INFO, rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+          line=__LINE__, &
+          file=__FILE__)) &
+          return  ! bail out
+        Time = esmf2fms_time(currTime-timeStep) ! i.e., startTime
+
+        call ESMF_LogWrite("MOM6 - doubling the coupling timestep", ESMF_LOGMSG_INFO, rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+          line=__LINE__, &
+          file=__FILE__)) &
+          return  ! bail out
+        Time_step_coupled = 2 * esmf2fms_time(timeStep)
+      endif
+    endif
+  endif
+
 
   !---------------
   ! Write diagnostics for import
@@ -1854,7 +1930,9 @@ subroutine ModelAdvance(gcomp, rc)
   !---------------
 
   if(profile_memory) call ESMF_VMLogMemInfo("Entering MOM update_ocean_model: ")
-  call update_ocean_model(Ice_ocean_boundary, ocean_state, ocean_public, Time, Time_step_coupled)
+  if (do_advance) then
+    call update_ocean_model(Ice_ocean_boundary, ocean_state, ocean_public, Time, Time_step_coupled)
+  endif
   if(profile_memory) call ESMF_VMLogMemInfo("Leaving MOM update_ocean_model: ")
 
   !---------------