add landice gridcomp

CHANGELOG.md

@@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Added landice grid comp
+
 ### Changed
 
 - More optimal distribution of tiles on processors for cubed-sphere tile space.

CMakeLists.txt

@@ -11,7 +11,7 @@ esma_add_library(${this}
   SRCS GEOS_LdasGridComp.F90
   SUBCOMPONENTS ${alldirs}
   SUBDIRS LDAS_Shared
-  DEPENDENCIES GEOSland_GridComp makebcs MAPL
+  DEPENDENCIES GEOSland_GridComp GEOSlandice_GridComp makebcs MAPL
   INCLUDES ${INC_ESMF})
 
 esma_add_subdirectory(GEOSldas_App)

GEOS_LdasGridComp.F90

@@ -14,6 +14,7 @@ module GEOS_LdasGridCompMod
   use GEOS_LandPertGridCompMod, only: LandPertSetServices => SetServices
   use GEOS_EnsGridCompMod, only: EnsSetServices => SetServices
   use GEOS_LandAssimGridCompMod, only: LandAssimSetServices => SetServices
+  use GEOS_LandiceGridCompMod,   only : LandiceSetServices  => SetServices
 
   use EASE_conv, only: ease_inverse
   use LDAS_TileCoordType, only: tile_coord_type , T_TILECOORD_STATE, TILECOORD_WRAP
@@ -50,6 +51,7 @@ module GEOS_LdasGridCompMod
 
   ! All children
   integer,allocatable :: LAND(:)
+  integer,allocatable :: LANDICE(:)
   integer,allocatable :: LANDPERT(:)
   integer,allocatable :: METFORCE(:)
   integer             :: ENSAVG, LANDASSIM
@@ -59,7 +61,7 @@ module GEOS_LdasGridCompMod
   logical :: land_assim
   logical :: mwRTM
   logical :: ensemble_forcing   ! switch between deterministic and ensemble forcing
-
+  logical :: with_landice
 contains
 
   !BOP
@@ -85,7 +87,7 @@ subroutine SetServices(gc, rc)
     character(len=ESMF_MAXSTR) :: Iam
     character(len=ESMF_MAXSTR) :: comp_name
     character(len=ESMF_MAXSTR) :: ensid_string,childname
-    character(len=ESMF_MAXSTR) :: LAND_ASSIM_STR, mwRTM_file, ENS_FORCING_STR
+    character(len=ESMF_MAXSTR) :: LAND_ASSIM_STR, mwRTM_file, ENS_FORCING_STR, WITH_LANDICE_STR
     integer                    :: ens_id_width
     ! Local variables
     type(T_TILECOORD_STATE), pointer :: tcinternal
@@ -94,7 +96,7 @@ subroutine SetServices(gc, rc)
     type(ESMF_Config) :: CF
     integer :: LSM_CHOICE
     integer :: FIRST_ENS_ID
-    
+
     ! Begin...
 
     ! Get my name and setup traceback handle
@@ -156,6 +158,12 @@ subroutine SetServices(gc, rc)
     VERIFY_(STATUS)
     land_assim = (trim(LAND_ASSIM_STR) /= 'NO')
 
+    call MAPL_GetResource ( MAPL, WITH_LANDICE_STR, Label="WITH_LANDICE:", DEFAULT="NO", RC=STATUS)
+    VERIFY_(STATUS)
+    WITH_LANDICE_STR =  ESMF_UtilStringUpperCase(WITH_LANDICE_STR, rc=STATUS)
+    VERIFY_(STATUS)
+    with_landice = (trim(WITH_LANDICE_STR) /= 'NO')
+
     call MAPL_GetResource ( MAPL, mwRTM_file, Label="LANDASSIM_INTERNAL_RESTART_FILE:", DEFAULT='', RC=STATUS)
     VERIFY_(STATUS)
     mwRTM = ( len_trim(mwRTM_file) /= 0 )
@@ -172,6 +180,7 @@ subroutine SetServices(gc, rc)
     endif
 
     allocate(LAND(NUM_ENSEMBLE),LANDPERT(NUM_ENSEMBLE))
+    if (with_landice) allocate(LANDICE(NUM_ENSEMBLE))
 
     ! ens_id_with = 2 + number of digits = total number of chars in ensid_string ("_eXXXX")
     !
@@ -212,6 +221,12 @@ subroutine SetServices(gc, rc)
        childname='LAND'//trim(ensid_string)
        LAND(i) = MAPL_AddChild(gc, name=childname, ss=LandSetServices, rc=status)
        VERIFY_(status)
+
+       if (with_landice) then
+          childname='LANDICE'//trim(ensid_string)
+          LANDICE(i) = MAPL_AddChild(gc, name=childname, ss=LandiceSetServices, rc=status)
+          VERIFY_(status)
+       endif
     enddo
 
     ENSAVG    = MAPL_AddChild(gc, name='ENSAVG', ss=EnsSetServices, rc=status)
@@ -224,19 +239,8 @@ subroutine SetServices(gc, rc)
 
     ! Connections
     do i=1,NUM_ENSEMBLE
-    ! -METFORCE-feeds-LANDPERT's-imports-
        k = 1
        if ( ensemble_forcing ) k = i
-       call MAPL_AddConnectivity(                                                  &
-            gc,                                                                    &
-            SHORT_NAME = ['Tair   ', 'Qair   ', 'Psurf  ', 'Rainf_C', 'Rainf  ',   &
-                          'Snowf  ', 'LWdown ', 'SWdown ', 'PARdrct', 'PARdffs',   &
-                          'Wind   ', 'RefH   '],                                   &
-            SRC_ID = METFORCE(k),                                                     &
-            DST_ID = LANDPERT(i),                                                  &
-            rc = status                                                            &
-            )
-       VERIFY_(status)
     ! -LANDPERT-feeds-LAND's-imports-
        call MAPL_AddConnectivity(                                                  &
             gc,                                                                    &
@@ -254,22 +258,7 @@ subroutine SetServices(gc, rc)
             rc = status                                                            &
             )
           VERIFY_(status)
-    ! -METFORCE-feeds-LAND's-imports-
-       call MAPL_AddConnectivity(                                                  &
-            gc,                                                                    &
-            SRC_NAME = ['Psurf', 'RefH ',                                          &
-                        'DUDP ', 'DUSV ', 'DUWT ', 'DUSD ', 'BCDP ', 'BCSV ',      &
-                        'BCWT ', 'BCSD ', 'OCDP ', 'OCSV ', 'OCWT ', 'OCSD ',      &
-                        'SUDP ', 'SUSV ', 'SUWT ', 'SUSD ', 'SSDP ', 'SSSV ' ],    &
-            SRC_ID = METFORCE(k),                                                     &
-            DST_NAME = ['PS  ', 'DZ  ',                                            &
-                        'DUDP', 'DUSV', 'DUWT', 'DUSD', 'BCDP', 'BCSV',            &
-                        'BCWT', 'BCSD', 'OCDP', 'OCSV', 'OCWT', 'OCSD',            &
-                        'SUDP', 'SUSV', 'SUWT', 'SUSD', 'SSDP', 'SSSV' ],          &
-            DST_ID = LAND(i),                                                      &
-            rc = status                                                            &
-            )
-       VERIFY_(status)
+
     ! -LAND-feeds-LANDPERT's-imports-
        call MAPL_AddConnectivity(                                                                  &
             gc,                                                                                    &
@@ -367,6 +356,8 @@ subroutine Initialize(gc, import, export, clock, rc)
     ! MAPL variables
     type(MAPL_LocStream) :: surf_locstream
     type(MAPL_LocStream) :: land_locstream
+    type(MAPL_LocStream) :: landice_locstream
+    type(MAPL_LocStream) :: force_locstream
     type(MAPL_MetaComp), pointer :: MAPL=>null() ! GC's MAPL obj
     type(MAPL_MetaComp), pointer :: CHILD_MAPL=>null() ! Child's MAPL obj
 
@@ -382,7 +373,7 @@ subroutine Initialize(gc, import, export, clock, rc)
     character(len=ESMF_MAXSTR) :: LAND_PARAMS 
     character(len=ESMF_MAXSTR) :: grid_type 
 
-    integer :: total_nt,land_nt_local,i,j
+    integer :: total_nt, land_nt_local, i, j
     real, pointer :: LandTileLats(:)
     real, pointer :: LandTileLons(:)
     integer, pointer :: local_id(:)
@@ -564,11 +555,32 @@ subroutine Initialize(gc, import, export, clock, rc)
     call MAPL_LocStreamCreate(                                                  &
          land_locstream,                                                        &
          surf_locstream,                                                        &
-         name=gcnames(LAND(1)),                                                    &
+         name=gcnames(LAND(1)),                                                 &
          mask=[MAPL_LAND],                                                      &
          rc=status                                                              &
          )
     VERIFY_(status)
+
+    if (with_landice) then
+       call MAPL_LocStreamCreate(                                                &
+          force_locstream,                                                       &
+          surf_locstream,                                                        &
+          name=gcnames(METFORCE(1)),                                             &
+          mask=[MAPL_LAND, MAPL_LANDICE],                                        &
+          rc=status                                                              &
+         )
+       VERIFY_(status)
+
+       call MAPL_LocStreamCreate(                                                &
+          landice_locstream,                                                     &
+          surf_locstream,                                                        &
+          name=gcnames(LANDICE(1)),                                              &
+          mask=[MAPL_LANDICE],                                                   &
+          rc=status                                                              &
+         )
+       VERIFY_(status)
+    endif
+
     call MAPL_TimerOff(MAPL, "-LocStreamCreate")
     ! Convert LAND's LocStream to LDAS' tile_coord and save it in the GridComp
     ! -get-tile-information-from-land's-locstream-
@@ -705,10 +717,16 @@ subroutine Initialize(gc, import, export, clock, rc)
     do i = 1, NUM_ENSEMBLE
        call MAPL_GetObjectFromGC(gcs(METFORCE(i)), CHILD_MAPL, rc=status)
        VERIFY_(status) ! CHILD = METFORCE
-       call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status)
-       VERIFY_(status)
+       if ( with_landice) then
+          call MAPL_Set(CHILD_MAPL, LocStream=force_locstream, rc=status)
+          VERIFY_(status)
+       else
+          call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status)
+          VERIFY_(status)
+       endif
        call ESMF_UserCompSetInternalState(gcs(METFORCE(i)), 'TILE_COORD', tcwrap, status)
        VERIFY_(status)
+
        ! exit after i=1 if using deterministic forcing
        if (.not. ensemble_forcing) exit
     enddo
@@ -723,15 +741,25 @@ subroutine Initialize(gc, import, export, clock, rc)
        VERIFY_(status)
        call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status)
        VERIFY_(status)
+
+       if (with_landice) then
+          call MAPL_GetObjectFromGC(gcs(LANDICE(i)), CHILD_MAPL, rc=status)
+          VERIFY_(status)
+          call MAPL_Set(CHILD_MAPL, LocStream=landice_locstream, rc=status)
+          VERIFY_(status)
+       endif
+
        call MAPL_GetObjectFromGC(gcs(LANDPERT(i)), CHILD_MAPL, rc=status)
        VERIFY_(status) ! CHILD = LANDPERT
        call MAPL_Set(CHILD_MAPL, LocStream=land_locstream, rc=status)
        VERIFY_(status)
+
        ! Add LAND's tile_coord to children's GridComps
        call ESMF_UserCompSetInternalState(gcs(LAND(i)), 'TILE_COORD', tcwrap, status)
        VERIFY_(status)
        call ESMF_UserCompSetInternalState(gcs(LANDPERT(i)), 'TILE_COORD', tcwrap, status)
        VERIFY_(status)
+
     enddo
 
     if (land_assim .or. mwRTM) then
@@ -822,7 +850,7 @@ subroutine Run(gc, import, export, clock, rc)
     type(MAPL_MetaComp), pointer :: MAPL
 
     ! Misc variables
-    integer :: igc,i, ens_id, FIRST_ENS_ID, ens_id_width
+    integer :: igc, i, ens_id, FIRST_ENS_ID, ens_id_width, k
     logical :: IAmRoot
     integer :: LSM_CHOICE
     type (ESMF_Field)                         :: field
@@ -893,13 +921,30 @@ subroutine Run(gc, import, export, clock, rc)
     do i = 1, NUM_ENSEMBLE
        igc = METFORCE(i)
        call MAPL_TimerOn(MAPL, gcnames(igc))
-       call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, userRC=status)
+       call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status)
        VERIFY_(status)
        call MAPL_TimerOff(MAPL, gcnames(igc))
        ! exit after i=1 if using deterministic forcing
        if (.not. ensemble_forcing) exit
     enddo
 
+    ! distribute force. ( export of focrce to the import of land, landpert and landice)
+    do i = 1, NUM_ENSEMBLE
+       k = 1
+       if (ensemble_forcing) k = i
+       igc = METFORCE(k)
+       call MAPL_TimerOn(MAPL, gcnames(igc))
+
+       call ESMF_GridCompRun(gcs(igc), importState=gex(igc), exportState=gim(LAND(i)),       clock=clock, phase=2, userRC=status)
+       VERIFY_(status)
+       call ESMF_GridCompRun(gcs(igc), importState=gex(igc), exportState=gim(LANDPERT(i)),   clock=clock, phase=3, userRC=status)
+       VERIFY_(status)
+       if (with_landice) then
+          call ESMF_GridCompRun(gcs(igc), importState=gex(igc), exportState=gim(LANDICE(i)),  clock=clock, phase=4, userRC=status)
+          VERIFY_(status)
+       endif
+       call MAPL_TimerOff(MAPL, gcnames(igc))
+    enddo
 
     do i  = 1,NUM_ENSEMBLE
 
@@ -927,6 +972,16 @@ subroutine Run(gc, import, export, clock, rc)
        VERIFY_(status)
        call MAPL_TimerOff(MAPL, gcnames(igc))
 
+       if (with_landice) then
+          igc = LANDICE(i)
+          call MAPL_TimerOn(MAPL, gcnames(igc))
+          call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=1, userRC=status)
+          VERIFY_(status)
+          call ESMF_GridCompRun(gcs(igc), importState=gim(igc), exportState=gex(igc), clock=clock, phase=2, userRC=status)
+          VERIFY_(status)
+          call MAPL_TimerOff(MAPL, gcnames(igc))
+       endif
+
        ! ApplyPrognPert - moved: now before calculating ensemble average that is picked up by land analysis and HISTORY; reichle 28 May 2020 
        igc = LANDPERT(i)
        call MAPL_TimerOn(MAPL, gcnames(igc))

GEOSldas_App/ldas_setup

@@ -110,7 +110,7 @@ class LDASsetup:
         self.out_path   = None
         self.inpdir   = None
         self.exefyl = None
-        self.islocal = False
+        self.isZoomIn = False
         self.catch = ''
         self.has_mwrtm = False
         self.has_vegopacity = False
@@ -740,18 +740,16 @@ class LDASsetup:
         tmp_f2g_file = tempfile.NamedTemporaryFile(delete=False) 
         cmd = self.bindir +'/preprocess_ldas.x c_f2g ' + tile + ' ' + self.domain_def.name + ' '+ self.out_path + ' ' + catchment_def + ' ' + exp_id + ' ' + _y4m2d2h2m2 + ' '+ dzsf + ' ' + tmp_f2g_file.name
 
-        print ('Creating f2g file: '+ tmp_f2g_file.name +'....\n')
+        print ('Creating f2g file if necessary: '+ tmp_f2g_file.name +'....\n')
         print ("cmd:   " + cmd)
         sp.call(shlex.split(cmd))
         # check if it is local or global
-        with open(tmp_f2g_file.name) as f2gfile :
-            head=[next(f2gfile) for x in range(2)]
-            if(head[0].strip() != head[1].strip()) :
-                self.islocal= True
+        if os.path.getsize(tmp_f2g_file.name) !=0 :
+           self.isZoomIn= True
         #os.remove(self.domain_def.name)
 
         # update tile domain
-        if self.islocal:
+        if self.isZoomIn:
             newlocalTile = tile+'.domain'
             print ("\nCreating local tile file :"+ newlocalTile)
             print ("\n by excluding land type MAPL_Land_ExcludeFromDomain=1100...\n")
@@ -780,7 +778,7 @@ class LDASsetup:
             bcstmp=bcstmp+[self.bcsdir+'/'+os.path.basename(bcf)]
         bcs=bcstmp
 
-        if self.islocal:
+        if self.isZoomIn:
            print ("Creating the boundary files for the simulation domain...\n")
            bcs_tmp=[]
            for bcf in bcs :
@@ -938,7 +936,7 @@ class LDASsetup:
             # catchment restart file
             if os.path.isfile(catchRstFile) :
                 catchLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2
-                if self.islocal :
+                if self.isZoomIn :
                     print( "Creating local catchment restart file... \n")
                     cmd=self.bindir +'/preprocess_ldas.x c_localcatchrst '+ catchRstFile +' ' + catchLocal + ' '+ tmp_f2g_file.name
                     print ("cmd:  "+cmd)
@@ -956,7 +954,7 @@ class LDASsetup:
             # vegdyn restart file
             if os.path.isfile(vegdynRstFile) :
                 vegdynLocal = self.rstdir+ensdir +'/'+self.rqdExeInp['EXP_ID']+'.vegdyn_internal_rst'
-                if self.islocal :
+                if self.isZoomIn :
                     print ("Creating the local veg restart file... \n")
                     cmd=self.bindir + '/preprocess_ldas.x c_localvegrst '+ vegdynRstFile +' ' + vegdynLocal + ' '+ tmp_f2g_file.name
                     print ("cmd:   " + cmd)
@@ -991,7 +989,7 @@ class LDASsetup:
         if self.has_mwrtm :
            mwRTMRstFile = self.mwrtm_file
            mwRTMLocal = self.bcsdir+'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.ldas_mwRTMparam.'+y4m2d2_h2m2+'z.nc4'
-           if self.islocal :
+           if self.isZoomIn :
               print ("Creating the local mwRTM restart file... \n")
               cmd= self.bindir +'/preprocess_ldas.x c_localmwrtmrst '+ mwRTMRstFile +' ' + mwRTMLocal + ' '+ tmp_f2g_file.name
 
@@ -1376,7 +1374,7 @@ class LDASsetup:
                   MY_EXPDOMAIN       = self.rqdExeInp['EXP_DOMAIN'],
                   MY_LOGFILE         = my_logfile,
                   MY_ERRFILE         = my_errfile,
-                  MY_MODEL           = self.catch,
+                  MY_LANDMODEL       = self.catch,
                   MY_POSTPROC_HIST   = str(self.rqdExeInp['POSTPROC_HIST']),
                   MY_FIRST_ENS_ID    = str(self.first_ens_id),
                   MY_LADAS_COUPLING  = str(self.ladas_coupling),