A new driver/coupler for GEOS-ESM

.gitignore

@@ -12,3 +12,6 @@ doc/build
 # Ignore testsuite file/directories
 testsuite*
 caselist*
+
+# icepack on its own repo
+icepack/

cicecore/cicedyn/dynamics/ice_dyn_evp.F90

@@ -418,13 +418,41 @@ subroutine evp (dt)
       !-----------------------------------------------------------------
 
       call stack_fields(tmass, aice_init, cdn_ocn, fld3)
+      call ice_timer_start(timer_bound)
       call ice_HaloUpdate (fld3,             halo_info, &
                            field_loc_center, field_type_scalar)
+      call ice_timer_stop(timer_bound)
+      call unstack_fields(fld3, tmass, aice_init, cdn_ocn)
+#ifdef GEOSCOUPLED
+      if (grid_ice == 'CD' .or. grid_ice == 'C') then
+         call stack_fields(uocn, ss_tltx, fld2)
+         call ice_timer_start(timer_bound)
+         call ice_HaloUpdate (fld2,           halo_info, &
+                              field_loc_Eface, field_type_vector)
+         call ice_timer_stop(timer_bound)
+         call unstack_fields(fld2, uocn, ss_tltx)
+         call stack_fields(vocn, ss_tlty, fld2)
+         call ice_timer_start(timer_bound)
+         call ice_HaloUpdate (fld2,           halo_info, &
+                              field_loc_Nface, field_type_vector)
+         call ice_timer_stop(timer_bound)
+         call unstack_fields(fld2, vocn, ss_tlty)
+      else
+         call stack_fields(uocn, vocn, ss_tltx, ss_tlty, fld4)
+         call ice_timer_start(timer_bound)
+         call ice_HaloUpdate (fld4,             halo_info, &
+                           field_loc_NEcorner, field_type_vector)
+         call ice_timer_stop(timer_bound)
+         call unstack_fields(fld4, uocn, vocn, ss_tltx, ss_tlty)
+      endif
+#else
       call stack_fields(uocn, vocn, ss_tltx, ss_tlty, fld4)
+      call ice_timer_start(timer_bound)
       call ice_HaloUpdate (fld4,             halo_info, &
                            field_loc_center, field_type_vector)
-      call unstack_fields(fld3, tmass, aice_init, cdn_ocn)
+      call ice_timer_stop(timer_bound)
       call unstack_fields(fld4, uocn, vocn, ss_tltx, ss_tlty)
+#endif
 
       call grid_average_X2Y('S', tmass    , 'T'          , umass   , 'U')
       call grid_average_X2Y('S', aice_init, 'T'          , aiU     , 'U')

cicecore/cicedyn/general/ice_flux.F90

@@ -184,6 +184,24 @@ module ice_flux
          fsensn_f   , & ! sensible heat flux (W m-2)
          flatn_f        ! latent heat flux (W m-2)
 
+
+#ifdef GEOSCOUPLED
+      ! in from atmosphere when coupling to GEOS
+      real (kind=dbl_kind), & 
+         dimension (:,:,:,:), allocatable, public :: &
+         evapn_f,      & ! evaporation/sublimation (kg m-2 s-1)
+         dflatndts_f,  & ! derivative of latent flux w.r.t. Ts
+         dfsurfndts_f    ! derivative of surface flux w.r.t. Ts
+
+      real (kind=dbl_kind), & 
+         dimension (:,:,:), allocatable, public :: &
+         druvr     , & !  vis uvr flux, direct (W m-2)
+         dfuvr     , & !  vis uvr flux, diffuse (W m-2)
+         drpar     , & !  vis par flux, direct (W m-2)
+         dfpar         !  vis par flux, diffuse (W m-2)
+#endif
+
+
        ! in from atmosphere
 
       real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
@@ -255,6 +273,10 @@ module ice_flux
          fhocn   , & ! net heat flux to ocean (W/m^2)
          fsloss  , & ! rate of snow loss to leads (kg/m^2/s)
          fswthru , & ! shortwave penetrating to ocean (W/m^2)
+         fswthru_uvrdr , & ! vis dir uvr SW penetrating to ocean (W/m^2)
+         fswthru_uvrdf , & ! vis dif uvr SW penetrating to ocean (W/m^2)
+         fswthru_pardr , & ! nir dir par SW penetrating to ocean (W/m^2)
+         fswthru_pardf , & ! nir dif par SW penetrating to ocean (W/m^2)
          fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2)
          fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2)
          fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2)
@@ -500,6 +522,10 @@ subroutine alloc_flux
          fhocn      (nx_block,ny_block,max_blocks), & ! net heat flux to ocean (W/m^2)
          fsloss     (nx_block,ny_block,max_blocks), & ! rate of snow loss to leads (kg/m^2/s)
          fswthru    (nx_block,ny_block,max_blocks), & ! shortwave penetrating to ocean (W/m^2)
+         fswthru_uvrdr (nx_block,ny_block,max_blocks), & ! vis dir uvr SW penetrating to ocean (W/m^2)
+         fswthru_uvrdf (nx_block,ny_block,max_blocks), & ! vis dir uvr SW penetrating to ocean (W/m^2)
+         fswthru_pardr (nx_block,ny_block,max_blocks), & ! vis dir par SW penetrating to ocean (W/m^2)
+         fswthru_pardf (nx_block,ny_block,max_blocks), & ! vis dir par SW penetrating to ocean (W/m^2)
          fswthru_vdr (nx_block,ny_block,max_blocks), & ! vis dir shortwave penetrating to ocean (W/m^2)
          fswthru_vdf (nx_block,ny_block,max_blocks), & ! vis dif shortwave penetrating to ocean (W/m^2)
          fswthru_idr (nx_block,ny_block,max_blocks), & ! nir dir shortwave penetrating to ocean (W/m^2)
@@ -568,6 +594,15 @@ subroutine alloc_flux
          fcondtopn_f(nx_block,ny_block,ncat,max_blocks), & ! downward cond flux at top surface (W m-2)
          fsensn_f   (nx_block,ny_block,ncat,max_blocks), & ! sensible heat flux (W m-2)
          flatn_f    (nx_block,ny_block,ncat,max_blocks), & ! latent heat flux (W m-2)
+#ifdef GEOSCOUPLED
+         evapn_f    (nx_block,ny_block,ncat,max_blocks), & ! evaporative water flux (kg/m^2/s) by atmosphere model
+         dflatndts_f (nx_block,ny_block,ncat,max_blocks), & ! derivative of flatn with respect to Ts 
+         dfsurfndts_f(nx_block,ny_block,ncat,max_blocks), & ! derivative of fsurfn with respect to Ts 
+         druvr      (nx_block,ny_block,max_blocks), & ! 
+         dfuvr      (nx_block,ny_block,max_blocks), & ! 
+         drpar      (nx_block,ny_block,max_blocks), & ! 
+         dfpar      (nx_block,ny_block,max_blocks), & ! 
+#endif
          meltsn     (nx_block,ny_block,ncat,max_blocks), & ! snow melt in category n (m)
          melttn     (nx_block,ny_block,ncat,max_blocks), & ! top melt in category n (m)
          meltbn     (nx_block,ny_block,ncat,max_blocks), & ! bottom melt in category n (m)
@@ -794,6 +829,10 @@ subroutine init_coupler_flux
       fswthru_vdf (:,:,:) = c0
       fswthru_idr (:,:,:) = c0
       fswthru_idf (:,:,:) = c0
+      fswthru_uvrdr (:,:,:) = c0
+      fswthru_uvrdf (:,:,:) = c0
+      fswthru_pardr (:,:,:) = c0
+      fswthru_pardf (:,:,:) = c0
       fresh_da(:,:,:) = c0    ! data assimilation
       fsalt_da(:,:,:) = c0
       flux_bio (:,:,:,:) = c0 ! bgc
@@ -855,6 +894,8 @@ subroutine init_flux_atm
       ! strairxT(:,:,:) = 0.15_dbl_kind
       ! strairyT(:,:,:) = 0.15_dbl_kind
 
+      fsurf   (:,:,:) = c0
+      fcondtop(:,:,:) = c0
       fsens   (:,:,:) = c0
       flat    (:,:,:) = c0
       fswabs  (:,:,:) = c0

cicecore/cicedyn/general/ice_init.F90

@@ -110,6 +110,9 @@ subroutine input_data
           bathymetry_file, use_bathymetry, &
           bathymetry_format, kmt_type, &
           grid_type, grid_format, &
+#ifdef GEOSCOUPLED
+          grid_subtype, & 
+#endif
           grid_ice, grid_ice_thrm, grid_ice_dynu, grid_ice_dynv, &
           grid_ocn, grid_ocn_thrm, grid_ocn_dynu, grid_ocn_dynv, &
           grid_atm, grid_atm_thrm, grid_atm_dynu, grid_atm_dynv, &
@@ -212,6 +215,9 @@ subroutine input_data
 
       namelist /grid_nml/ &
         grid_format,    grid_type,       grid_file,     kmt_file,       &
+#ifdef GEOSCOUPLED
+        grid_subtype,                                                   &
+#endif
         bathymetry_file, use_bathymetry, nfsd,          bathymetry_format, &
         ncat,           nilyr,           nslyr,         nblyr,          &
         kcatbound,      gridcpl_file,    dxrect,        dyrect,         &
@@ -377,6 +383,9 @@ subroutine input_data
       ice_ic       = 'default'    ! latitude and sst-dependent
       grid_format  = 'bin'        ! file format ('bin'=binary or 'nc'=netcdf)
       grid_type    = 'rectangular'! define rectangular grid internally
+#ifdef GEOSCOUPLED
+      grid_subtype  = 'other'  ! not belonging to a specific subtype 
+#endif
       grid_file    = 'unknown_grid_file'
       grid_ice     = 'B'          ! underlying grid system
       grid_atm     = 'A'          ! underlying atm forcing/coupling grid
@@ -985,6 +994,9 @@ subroutine input_data
       call broadcast_scalar(latrefrect,           master_task)
       call broadcast_scalar(close_boundaries,     master_task)
       call broadcast_scalar(grid_type,            master_task)
+#ifdef GEOSCOUPLED
+      call broadcast_scalar(grid_subtype,         master_task)
+#endif
       call broadcast_scalar(grid_ice,             master_task)
       call broadcast_scalar(grid_ocn,             master_task)
       call broadcast_scalar(grid_atm,             master_task)

cicecore/cicedyn/general/ice_state.F90

@@ -133,6 +133,10 @@ module ice_state
          vicen_init  , & ! initial ice volume (m), for linear ITD
          vsnon_init      ! initial snow volume (m), for aerosol
 
+      real (kind=dbl_kind), &
+         dimension(:,:,:,:), allocatable, public :: &
+         Tsfcn_init ! initial ice surface temperature (degC)   
+
 !=======================================================================
 
       contains
@@ -172,7 +176,8 @@ subroutine alloc_state
          vsnon     (nx_block,ny_block,ncat,max_blocks) , & ! volume per unit area of snow (m)
          aicen_init(nx_block,ny_block,ncat,max_blocks) , & ! initial ice concentration, for linear ITD
          vicen_init(nx_block,ny_block,ncat,max_blocks) , & ! initial ice volume (m), for linear ITD
-         vsnon_init(nx_block,ny_block,ncat,max_blocks) , & ! initial snow volume (m), for aerosol
+         vsnon_init(nx_block,ny_block,ncat,max_blocks) , & ! initial snow volume (m), for aerosol 
+         Tsfcn_init(nx_block,ny_block,ncat,max_blocks) , & ! initial snow/ice surface temperature(degC)
          trcr      (nx_block,ny_block,ntrcr,max_blocks) , & ! ice tracers: 1: surface temperature of ice/snow (C)
          trcrn     (nx_block,ny_block,ntrcr,ncat,max_blocks) , & ! tracers: 1: surface temperature of ice/snow (C)
          stat=ierr)

cicecore/cicedyn/general/ice_step_mod.F90

@@ -60,7 +60,7 @@ subroutine save_init
 ! saves initial values for aice, aicen, vicen, vsnon
 
       use ice_state, only: aice, aicen, aice_init, aicen_init, &
-          vicen, vicen_init, vsnon, vsnon_init
+          vicen, vicen_init, vsnon, vsnon_init, trcrn, Tsfcn_init
 
       character(len=*), parameter :: subname = '(save_init)'
 
@@ -74,6 +74,7 @@ subroutine save_init
          aicen_init = aicen
          vicen_init = vicen
          vsnon_init = vsnon
+         Tsfcn_init = trcrn(:,:,1,:,:)
 
       end subroutine save_init
 
@@ -222,6 +223,7 @@ subroutine step_therm1 (dt, iblk)
           Cdn_atm, Cdn_atm_skin, Cdn_atm_floe, Cdn_atm_rdg, Cdn_atm_pond, &
           hfreebd, hdraft, hridge, distrdg, hkeel, dkeel, lfloe, dfloe, &
           fswsfcn, fswintn, Sswabsn, Iswabsn, meltsliqn, meltsliq, &
+          fswthrun_uvrdr, fswthrun_uvrdf, fswthrun_pardr, fswthrun_pardf, &
           fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf
       use ice_calendar, only: yday
       use ice_domain_size, only: ncat, nilyr, nslyr, n_iso, n_aero
@@ -232,12 +234,20 @@ subroutine step_therm1 (dt, iblk)
           frain, Tair, strairxT, strairyT, fsurf, fcondtop, fsens, &
           flat, fswabs, flwout, evap, evaps, evapi, Tref, Qref, Uref, fresh, fsalt, fhocn, &
           fswthru, fswthru_vdr, fswthru_vdf, fswthru_idr, fswthru_idf, &
+          fswthru_uvrdr, fswthru_uvrdf, fswthru_pardr, fswthru_pardf, &
           meltt, melts, meltb, congel, snoice, &
+#ifdef GEOSCOUPLED
+          dfsurfndts_f, dflatndts_f,  &
+#endif
           flatn_f, fsensn_f, fsurfn_f, fcondtopn_f, &
           send_i2x_per_cat, fswthrun_ai, dsnow
       use ice_flux_bgc, only: dsnown, faero_atm, faero_ocn, fiso_atm, fiso_ocn, &
           Qa_iso, Qref_iso, fiso_evap, HDO_ocn, H2_16O_ocn, H2_18O_ocn
       use ice_grid, only: lmask_n, lmask_s, tmask
+#ifdef GEOSCOUPLED
+      use ice_communicate, only: my_task
+      use ice_grid, only: opmask
+#endif
       use ice_state, only: aice, aicen, aicen_init, vicen_init, &
           vice, vicen, vsno, vsnon, trcrn, vsnon_init
 #ifdef CICE_IN_NEMO
@@ -387,7 +397,11 @@ subroutine step_therm1 (dt, iblk)
             enddo
          endif ! tr_aero
 
+#ifdef GEOSCOUPLED
+         if (tmask(i,j,iblk) .or. opmask(i,j,iblk)) then
+#else
          if (tmask(i,j,iblk)) then
+#endif
 
          call icepack_step_therm1(dt=dt, ncat=ncat,            &
                       nilyr=nilyr, nslyr=nslyr,                &
@@ -485,6 +499,10 @@ subroutine step_therm1 (dt, iblk)
                       fswthrun_vdf = fswthrun_vdf (i,j,:,iblk),&
                       fswthrun_idr = fswthrun_idr (i,j,:,iblk),&
                       fswthrun_idf = fswthrun_idf (i,j,:,iblk),&
+                      fswthrun_uvrdr = fswthrun_uvrdr (i,j,:,iblk),&
+                      fswthrun_uvrdf = fswthrun_uvrdf (i,j,:,iblk),&
+                      fswthrun_pardr = fswthrun_pardr (i,j,:,iblk),&
+                      fswthrun_pardf = fswthrun_pardf (i,j,:,iblk),&
                       fswabs       = fswabs      (i,j,  iblk), &
                       flwout       = flwout      (i,j,  iblk), &
                       Sswabsn      = Sswabsn     (i,j,:,:,iblk), &
@@ -505,6 +523,10 @@ subroutine step_therm1 (dt, iblk)
                       fswthru_vdf  = fswthru_vdf (i,j,  iblk), &
                       fswthru_idr  = fswthru_idr (i,j,  iblk), &
                       fswthru_idf  = fswthru_idf (i,j,  iblk), &
+                      fswthru_uvrdr  = fswthru_uvrdr (i,j,  iblk), &
+                      fswthru_uvrdf  = fswthru_uvrdf (i,j,  iblk), &
+                      fswthru_pardr  = fswthru_pardr (i,j,  iblk), &
+                      fswthru_pardf  = fswthru_pardf (i,j,  iblk), &
                       flatn_f      = flatn_f     (i,j,:,iblk), &
                       fsensn_f     = fsensn_f    (i,j,:,iblk), &
                       fsurfn_f     = fsurfn_f    (i,j,:,iblk), &
@@ -537,7 +559,13 @@ subroutine step_therm1 (dt, iblk)
                       lmask_s      = lmask_s     (i,j,  iblk), &
                       mlt_onset    = mlt_onset   (i,j,  iblk), &
                       frz_onset    = frz_onset   (i,j,  iblk), &
-                      yday=yday, prescribed_ice=prescribed_ice)
+#ifdef GEOSCOUPLED
+                      dfsurfdt     = dfsurfndts_f(i,j,:,iblk), &
+                      dflatdt      = dflatndts_f (i,j,:,iblk), &
+#endif
+                      yday         = yday,                     &
+                      prescribed_ice=prescribed_ice            ) 
+
 
       !-----------------------------------------------------------------
       ! handle per-category i2x fields, no merging
@@ -622,6 +650,9 @@ subroutine step_therm2 (dt, iblk)
       use ice_flux_bgc, only: flux_bio, faero_ocn, &
           fiso_ocn, HDO_ocn, H2_16O_ocn, H2_18O_ocn
       use ice_grid, only: tmask
+#ifdef GEOSCOUPLED
+      use ice_grid, only: opmask
+#endif
       use ice_state, only: aice, aicen, aice0, trcr_depend, &
           aicen_init, vicen_init, trcrn, vicen, vsnon, &
           trcr_base, n_trcr_strata, nt_strata
@@ -673,7 +704,11 @@ subroutine step_therm2 (dt, iblk)
       do j = jlo, jhi
       do i = ilo, ihi
 
+#ifdef GEOSCOUPLED
+         if (tmask(i,j,iblk) .or. opmask(i,j,iblk)) then
+#else
          if (tmask(i,j,iblk)) then
+#endif
 
          ! significant wave height for FSD
          if (tr_fsd) &
@@ -1273,13 +1308,18 @@ subroutine step_radiation (dt, iblk)
       use ice_arrays_column, only: ffracn, dhsn, &
           fswsfcn, fswintn, fswpenln, Sswabsn, Iswabsn, &
           fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf, &
+          fswthrun_uvrdr, fswthrun_uvrdf, fswthrun_pardr, fswthrun_pardf,   &
           albicen, albsnon, albpndn, &
           alvdrn, alidrn, alvdfn, alidfn, apeffn, trcrn_sw, snowfracn, &
           swgrid, igrid
       use ice_calendar, only: calendar_type, days_per_year, nextsw_cday, yday, msec
       use ice_domain_size, only: ncat, n_aero, nilyr, nslyr, n_zaero, n_algae, nblyr
       use ice_flux, only: swvdr, swvdf, swidr, swidf, coszen, fsnow
       use ice_grid, only: TLAT, TLON, tmask
+#ifdef GEOSCOUPLED
+      use ice_grid, only: opmask
+      use ice_flux, only: druvr, dfuvr, drpar, dfpar 
+#endif
       use ice_state, only: aicen, vicen, vsnon, trcrn
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_sw
       use ice_communicate, only: my_task
@@ -1379,8 +1419,11 @@ subroutine step_radiation (dt, iblk)
                enddo
             endif
          enddo
-
+#ifdef GEOSCOUPLED
+         if (tmask(i,j,iblk) .or. opmask(i,j,iblk)) then
+#else
          if (tmask(i,j,iblk)) then
+#endif
 
             call icepack_step_radiation (dt=dt,                               &
                          swgrid=swgrid(:),        igrid=igrid(:),             &
@@ -1404,6 +1447,10 @@ subroutine step_radiation (dt, iblk)
                          sec=msec,                                             &
                          swvdr    =swvdr    (i,j    ,iblk), swvdf   =swvdf   (i,j    ,iblk), &
                          swidr    =swidr    (i,j    ,iblk), swidf   =swidf   (i,j    ,iblk), &
+#ifdef GEOSCOUPLED
+                         druvr    =druvr    (i,j    ,iblk), dfuvr   =dfuvr   (i,j    ,iblk), &
+                         drpar    =drpar    (i,j    ,iblk), dfpar   =dfpar   (i,j    ,iblk), &
+#endif
                          coszen   =coszen   (i,j    ,iblk), fsnow   =fsnow   (i,j    ,iblk), &
                          alvdrn   =alvdrn   (i,j,:  ,iblk), alvdfn  =alvdfn  (i,j,:  ,iblk), &
                          alidrn   =alidrn   (i,j,:  ,iblk), alidfn  =alidfn  (i,j,:  ,iblk), &
@@ -1413,6 +1460,10 @@ subroutine step_radiation (dt, iblk)
                          fswthrun_vdf =fswthrun_vdf (i,j,:  ,iblk), &
                          fswthrun_idr =fswthrun_idr (i,j,:  ,iblk), &
                          fswthrun_idf =fswthrun_idf (i,j,:  ,iblk), &
+                         fswthrun_uvrdr =fswthrun_uvrdr (i,j,:  ,iblk), &
+                         fswthrun_uvrdf =fswthrun_uvrdf (i,j,:  ,iblk), &
+                         fswthrun_pardr =fswthrun_pardr (i,j,:  ,iblk), &
+                         fswthrun_pardf =fswthrun_pardf (i,j,:  ,iblk), &
                          fswpenln=fswpenln(i,j,:,:,iblk), &
                          Sswabsn  =Sswabsn  (i,j,:,:,iblk), Iswabsn =Iswabsn (i,j,:,:,iblk), &
                          albicen  =albicen  (i,j,:  ,iblk), albsnon =albsnon (i,j,:  ,iblk), &
@@ -1444,6 +1495,7 @@ subroutine step_radiation (dt, iblk)
 
       end subroutine step_radiation
 
+
 !=======================================================================
 ! Ocean mixed layer calculation (internal to sea ice model).
 ! Allows heat storage in ocean for uncoupled runs.

cicecore/cicedyn/infrastructure/ice_blocks.F90

@@ -323,6 +323,12 @@ subroutine create_blocks(nx_global, ny_global, ew_boundary_type, &
          all_blocks(n)%jblock,   &
          all_blocks(n)%tripole
       enddo
+      write(nu_diag,'(*(i3, 1x))') ( i_global(all_blocks(i)%ihi,i) -  & 
+                                     i_global(all_blocks(i)%ilo,i)+1, & 
+                                     i=1,nblocks_x ) 
+      write(nu_diag,'(*(i3, 1x))') ( j_global(all_blocks((j-1)*nblocks_x+1)%jhi,(j-1)*nblocks_x+1) -  & 
+                                     j_global(all_blocks((j-1)*nblocks_x+1)%jlo,(j-1)*nblocks_x+1)+1, & 
+                                     j=1,nblocks_y ) 
       endif
    endif