Merge pull request #753 from GEOS-ESM/develop

Sync develop into main

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSlake_GridComp/GEOS_LakeGridComp.F90

@@ -1,5 +1,17 @@
 !  $Id$
 
+#define _DEALLOC(A) \
+    if(associated(A))then; \
+          if(MAPL_ShmInitialized)then; \
+              call MAPL_SyncSharedMemory(rc=STATUS); \
+              call MAPL_DeAllocNodeArray(A,rc=STATUS); \
+           else; \
+              deallocate(A,stat=STATUS); \
+           endif; \
+       _VERIFY(STATUS); \
+       NULLIFY(A); \
+    endif
+
 #include "MAPL_Generic.h"
 
 
@@ -25,7 +37,13 @@ module GEOS_LakeGridCompMod
   integer, parameter :: NUM_SUBTILES = 2
 
   type lake_state
-       integer:: CHOOSEMOSFC
+     private
+     integer:: CHOOSEMOSFC
+     logical :: InitDone
+     logical, pointer :: mask(:) => null()
+     real :: tol_frice
+     character(len=ESMF_MAXSTR) :: sstfile
+     character(len=ESMF_MAXSTR) :: DataFrtFile
   end type lake_state
 
   type lake_state_wrap
@@ -854,6 +872,7 @@ subroutine SetServices ( GC, RC )
 
     allocate(mystate,stat=status)
     VERIFY_(status)
+    mystate%InitDone = .false.
     call MAPL_GetResource (MAPL, SURFRC, label = 'SURFRC:', default = 'GEOS_SurfaceGridComp.rc', RC=STATUS) ; VERIFY_(STATUS)
     SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS)
     call ESMF_ConfigLoadFile     (SCF,SURFRC,rc=status) ; VERIFY_(STATUS)
@@ -914,8 +933,6 @@ subroutine RUN1 ( GC, IMPORT, EXPORT, CLOCK, RC )
 
   type (MAPL_MetaComp), pointer   :: MAPL
   type (ESMF_State       )            :: INTERNAL
-  type (ESMF_Alarm       )            :: ALARM
-  type (ESMF_Config      )            :: CF
 
 ! pointers to export
 
@@ -1553,6 +1570,13 @@ subroutine LAKECORE(NT,RC)
    real, parameter :: LAKECAP      = (MAPL_RHOWTR*MAPL_CAPWTR*LAKEDEPTH    )
    real, parameter :: LAKEICECAP   = (MAPL_RHOWTR*MAPL_CAPWTR*LAKEICEDEPTH )
 
+   logical :: datalake
+   integer :: dlk
+   real, allocatable :: DATA_SST(:), DATA_FR(:)
+   character(len=ESMF_MAXSTR) :: maskfile
+   real, parameter :: Tfreeze = MAPL_TICE - 1.8
+   type(lake_state_wrap) :: wrap
+   type(lake_state), pointer :: mystate
 
 !  Begin...
 !----------
@@ -1713,6 +1737,80 @@ subroutine LAKECORE(NT,RC)
     if(associated(QST    )) QST     = 0.0
     if(associated(HLWUP  )) HLWUP   = ALW 
     if(associated(LWNDSRF)) LWNDSRF = LWDNSRF - ALW
+    ! datalake addition
+    !==================
+
+    ! check resource if datalake is enabled
+    call MAPL_GetResource ( MAPL, DLK, Label="DATALAKE:", DEFAULT=0, RC=STATUS)
+    VERIFY_(STATUS)
+    DATALAKE = (DLK /= 0)
+
+    call ESMF_UserCompGetInternalState(gc,'lake_private',wrap,status)
+    VERIFY_(status)
+    mystate => wrap%ptr
+
+    ! Initialize a mask where we could apply the SST/FR from Reynolds/Ostia
+    if (.not. associated(mystate%mask)) then
+       allocate(mystate%mask(NT), stat=status); VERIFY_(STATUS)
+       mystate%mask = .false.
+    end if
+
+    if (datalake) then
+
+        ! next section is done only once. We do it here since the Initalize
+        ! method of this component defaults to MAPL_GenericInitialize
+
+        if (.not. mystate%InitDone) then
+           mystate%InitDone = .true.
+           call MAPL_GetResource ( MAPL, mystate%sstfile, &
+                Label="LAKE_SST_FILE:", DEFAULT="sst.data", RC=STATUS)
+           VERIFY_(STATUS)
+           call MAPL_GetResource ( MAPL, mystate%dataFrtFile, &
+                Label="LAKE_FRT_FILE:", DEFAULT="fraci.data", RC=STATUS)
+           VERIFY_(STATUS)
+           call MAPL_GetResource ( MAPL, mystate%tol_frice, &
+                Label="LAKE_TOL_FRICE:", DEFAULT=1.0e-2, RC=STATUS)
+           VERIFY_(STATUS)
+
+           call MAPL_GetResource ( MAPL, MASKFILE, &
+                Label="DATALAKE_MASK_FILE:", DEFAULT="DataLakeMask.data", RC=STATUS)
+           VERIFY_(STATUS)
+
+           !ALT: For now we are reading the entire mask from a file
+           ! we might decide later to use a different strategy (for example
+           ! Santha suggested lat-lon based description)
+
+           call DataLakeReadMask(MAPL, mystate%mask, &
+                maskfile, RC=status)
+           VERIFY_(STATUS)
+        end if
+
+
+       allocate(DATA_SST(NT), DATA_FR(NT), stat=status); VERIFY_(STATUS)
+       call MAPL_ReadForcing(MAPL, 'SST', mystate%sstfile,&
+            CURRENT_TIME, DATA_SST, RC=STATUS)
+       VERIFY_(STATUS)
+       call MAPL_ReadForcing(MAPL, 'FR', mystate%dataFrtFile,&
+            CURRENT_TIME, DATA_FR, RC=STATUS)
+       VERIFY_(STATUS)
+
+       do I=1,NT
+          if(mystate%mask(i)) then
+             ! we are operating over the 'observed' lake: Great Lakes and Caspian Sea (set by mask)
+             TS(I,WATER) = DATA_SST(I)
+             TS(I,ICE)   = Tfreeze
+             if (data_fr(i) > mystate%tol_frice) then !have lake ice
+                FR(I,WATER) = 1.0-DATA_FR(I)
+                FR(I,ICE)   = DATA_FR(I)
+             else
+                FR(I,WATER) = 1.0
+                FR(I,ICE)   = 0.0
+             end if
+          end if
+       end do
+
+       deallocate(DATA_SST, DATA_FR)
+    end if
 
     do N=1,NUM_SUBTILES
        CFT = (CH(:,N)/CTATM)
@@ -1749,7 +1847,7 @@ subroutine LAKECORE(NT,RC)
 ! Update surface temperature and moisture
 !----------------------------------------
 
-       TS(:,N) = TS(:,N) + DTS
+       where(.not.mystate%mask) TS(:,N) = TS(:,N) + DTS
        DQS     = GEOS_QSAT(TS(:,N), PS, RAMP=0.0, PASCALS=.TRUE.) - QS(:,N)
        QS(:,N) = QS(:,N) + DQS  
 
@@ -1773,6 +1871,10 @@ subroutine LAKECORE(NT,RC)
 ! Update Ice fraction
 !--------------------
     do I=1,NT
+!$       if(mystate%mask(i)) cycle 
+       if(mystate%mask(i)) then
+          cycle 
+       end if
        if    (TS(I,ICE)>MAPL_TICE .and. FR(I,ICE)>0.0) then
           ! MELT
           FR(I,WATER) = 1.0
@@ -1919,6 +2021,54 @@ end subroutine ALBLAKEICE
 
 end subroutine RUN2
 
+subroutine DataLakeReadMask(mapl, msk, maskfile, rc)
+  implicit none
+  ! arguments
+  type (MAPL_MetaComp), pointer  :: MAPL
+  logical, intent(INOUT)           :: msk(:)
+  character(len=*), intent(IN)   :: maskfile
+  integer, optional, intent(OUT) :: rc
+
+  ! errlog vars
+  integer :: status
+  character(len=ESMF_MAXSTR), parameter :: Iam='DataLakeReadMask'
+
+  ! local vars
+  integer                       :: unit
+  integer, pointer              :: tilemask(:) => null()
+  type(ESMF_Grid)               :: TILEGRID
+  type(MAPL_LocStream)          :: LOCSTREAM
+  integer :: NT
+  real, allocatable :: imask(:)
+
+  ! this is the first attempt to read the mask. For now we support binary
+  call MAPL_Get(MAPL, LocStream=LOCSTREAM, RC=STATUS)
+  VERIFY_(STATUS)
+  call MAPL_LocStreamGet(LOCSTREAM, TILEGRID=TILEGRID, RC=STATUS)
+  VERIFY_(STATUS)
+
+  call MAPL_TileMaskGet(tilegrid, tilemask, rc=status)
+  VERIFY_(STATUS)
+
+  nt = size(msk)
+  allocate(imask(nt), stat=status)
+
+  unit = GETFILE( maskfile, form="unformatted", RC=STATUS )
+  VERIFY_(STATUS)
+
+  call MAPL_VarRead(unit, tilegrid, imask,  mask=tilemask, rc=status)
+  VERIFY_(STATUS)
+
+  call FREE_FILE(unit, RC=STATUS)
+  VERIFY_(STATUS)
+
+  msk = (imask /= 0.0)
+  deallocate(imask)
+  _DEALLOC(tilemask)
+
+  RETURN_(ESMF_SUCCESS)
+end subroutine DataLakeReadMask
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 end module GEOS_LakeGridCompMod

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSlake_GridComp/maskLakes.F90

@@ -0,0 +1,102 @@
+#define _ASSERT(cond) if(.not. cond) then; print*, 'ERROR at ', __LINE__; stop; endif
+program maskLakes
+  implicit none
+
+  integer, parameter :: HDR_SIZE=14
+  integer :: header(HDR_SIZE)
+  real :: hdr(HDR_SIZE)
+
+  integer :: nt, nl
+  integer :: ng
+  character(len=128), parameter :: filename='mask.bin'
+  character(len=128), parameter :: tilefile='tile.bin'
+  real, allocatable :: lats(:), lons(:)
+  integer, allocatable :: t(:)
+  real, allocatable :: buffer(:)
+  real, allocatable :: omask(:,:) ! mask on OSTIA grid
+  integer :: im, jm
+  integer :: i, j, n
+  integer :: status
+  integer :: unit
+  real, allocatable :: mask(:)
+  real :: x0, y0, dx, dy
+
+! read OSTIA mask
+  unit=10
+  open(unit=unit, file=filename, form='unformatted')
+  read(unit) hdr
+  header = nint(hdr)
+  im = header(13)
+  jm = header(14)
+  print *,'DEBUG: im/jm',im,jm
+
+  allocate(omask(im,jm), stat=status)
+  _ASSERT(status == 0)
+  read(unit) omask
+  close(unit)
+
+! process tile file
+  unit=11
+  open(unit=unit, file=tilefile, form='unformatted')
+  read(unit) nt
+  print *,'DEBUG ntiles ',nt
+  read(unit) ng
+  print *,'DEBUG ngrids ',ng
+  do i=1,ng
+     ! 3 blank read for im, jm, gridname
+     read(unit)
+     read(unit)
+     read(unit)
+  end do
+  allocate(buffer(nt), stat=status)
+  _ASSERT(status == 0)
+! get typetype
+  read(unit) buffer
+  t = nint(buffer)
+  nl = count(t==19) ! number of lake tiles
+  print *,'DEBUG lake points ',nl
+  allocate(lons(nl), lats(nl), mask(nl), stat=status)
+  _ASSERT(status == 0)
+  mask = 0.0
+
+! get X
+  read(unit) buffer
+  lons = pack(buffer, t==19)
+! get Y
+  read(unit) buffer
+  lats = pack(buffer, t==19)
+  close(unit)
+
+  ! the OSTIA grid origin is at the pole edge, dateline edge
+  x0 = -180.0
+  y0 = -90.0
+  dx = 360. / im
+  dy = 180. / jm
+
+!  convert tile lat/lon to Ostia grid indices
+  do n=1,nl
+     i = nint((lons(n)-x0)/dx - 0.5)
+     j = nint((lats(n)-y0)/dy -0.5)+1
+     i = mod(i+im,im)+1
+     if (i<=0 .or. j<=0) print *,i,j,n,lats(n),lons(n)
+     if (i>im .or. j>jm) print *,i,j,n,lats(n),lons(n)
+     _ASSERT(i>0 .and. i<im)
+     _ASSERT(j>0 .and. j<jm)
+
+     if(omask(i,j) /= 0.0) then 
+        mask(n) = 1.0
+        print *,'DEBUG: masked lake at ',n,lons(n),lats(n)
+     end if
+  end do
+
+! some sanity check/prints
+  print *,'done, processed ',nl,' lake points, found ',count(mask/=0.0)
+
+! ultimately, write mask
+  unit=12
+  open(unit=unit, file='lakemask.bin', form='unformatted')
+  write(unit) mask
+  close(unit)
+
+! all done
+end program maskLakes

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90

@@ -5003,7 +5003,7 @@ subroutine Driver ( RC )
 
         type(ESMF_Config)           :: CF
         type(MAPL_SunOrbit)         :: ORBIT
-        type(ESMF_Time)             :: CURRENT_TIME, StopTime, NextTime
+        type(ESMF_Time)             :: CURRENT_TIME, StopTime, NextTime, NextRecordTime
         type(ESMF_Time)             :: BEFORE
         type(ESMF_Time)             :: NOW
         type(ESMF_Time)             :: MODELSTART
@@ -5167,6 +5167,10 @@ subroutine Driver ( RC )
     ! Variables for FPAR
         ! --------------------------
     real   , allocatable, dimension (:,:)      :: parzone
+
+    logical           :: record
+    type(ESMF_Alarm)  :: RecordAlarm 
+
     character(len=ESMF_MAXSTR) :: Co2_CycleFile
 
         IAm=trim(COMP_NAME)//"::RUN2::Driver"
@@ -6843,7 +6847,14 @@ subroutine Driver ( RC )
 
        ! copy CN_restart vars to catch_internal_rst gkw: only do if stopping
        ! ------------------------------------------
-       if(NextTime == StopTime) then
+       record = .false.
+       call ESMF_ClockGetAlarm ( CLOCK, alarmname="RecordAlarm001", ALARM=RecordAlarm, RC=STATUS )
+       if (status == 0) then
+           call ESMF_AlarmGet( RecordAlarm, RingTime=NextRecordTime, _RC)
+           if (NextTime == NextRecordTime) record = .true.
+       endif
+
+       if(NextTime == StopTime  .or. record ) then
 
           call CN_exit(ntiles,nveg,nzone,ityp,fveg,cncol,var_col,cnpft,var_pft)     
           i = 1