379 continue integrating the ground water hydrology work into cable

src/offline/cable_checks.F90

@@ -84,7 +84,7 @@ MODULE cable_checks_module
       ESoil = [-0.0015, 0.0015], &
       TVeg = [-0.0003, 0.0003], &
       ECanop = [-0.0003, 0.0003], &
-      PotEvap = [-0.005, 0.005], &  !note should encompass Evap
+      PotEvap = [-0.005, 0.005], &  !note should encompass Evap  ! I have PotEvap = [-0.0006, 0.0006] if it makes a difference - rk4417
       ACond = [0.0, 1.0], &
       SoilWet = [-0.4, 1.2], &
       Albedo = [0.0, 1.0], &
@@ -133,12 +133,15 @@ MODULE cable_checks_module
       css = [700.0, 2200.0], &
       rhosoil = [300.0, 3000.0], &
       hyds = [5.0E-7, 8.5E-3], & ! vh_js ! sep14
+! MMY 8.5E-3->8.5 since hyds uses m/s, but hyds_vec uses mm/s
       rs20 = [0.0, 10.0], &
       sand = [0.0, 1.0], &
       sfc = [0.1, 0.5], &
       silt = [0.0, 1.0], &
       ssat = [0.35, 0.5], &
-      sucs = [-0.8, -0.03], &
+! 2 lines below changed by rk4417 - phase2
+!      sucs = [-0.8, -0.03],              & ! MMY@23Apr2023 keep this line commented since it works for CABLE-non GW
+      sucs = [30., 800.],                & ! MMY the range [-0.8, -0.03] doesn't suit for Mark Decker's version
       swilt = [0.05, 0.4], &
       froot = [0.0, 1.0], &
       zse = [0.0, 5.0], &
@@ -486,6 +489,8 @@ SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met,
 
     ! Local variables
     REAL(r_2), DIMENSION(:,:,:),POINTER, SAVE :: bwb         ! volumetric soil moisture
+! line below inserted by rk4417 - phase2
+    REAL(r_2), DIMENSION(:,:),POINTER, SAVE   :: bwb_gw ! volumetric gw soil moisture ! MMY
     REAL(r_2), DIMENSION(mp)                  :: delwb       ! change in soilmoisture
     ! b/w tsteps
     REAL, DIMENSION(mp)                  :: canopy_wbal !canopy water balance
@@ -494,12 +499,39 @@ SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met,
 
     IF(ktau==1) THEN
        ALLOCATE( bwb(mp,ms,2) )
-       ! initial vlaue of soil moisture
+       ALLOCATE( bwb_gw(mp,2) ) ! MMY  ! inserted by rk4417 - phase2
+       ! initial value of soil moisture
        bwb(:,:,1)=ssnow%wb
        bwb(:,:,2)=ssnow%wb
+       bwb_gw(:,1)=ssnow%GWwb ! MMY   ! 2 lines inserted by rk4417 - phase2
+       bwb_gw(:,2)=ssnow%GWwb ! MMY
        delwb(:) = 0.
     ELSE
        ! Calculate change in soil moisture b/w timesteps:
+
+      ! ________________ MMY, Water Balance Equation for GW_ON _______________
+      ! MMY to fix water balance bug when gw-on
+      IF ( cable_user%GW_MODEL) THEN ! MMY
+
+       IF(MOD(REAL(ktau),2.0)==1.0) THEN         ! if odd timestep
+          bwb(:,:,1)=ssnow%wb
+          bwb_gw(:,1)=ssnow%GWwb
+          DO k=1,mp           ! current smoist - prev tstep smoist
+             delwb(k) = ( SUM( (bwb(k,:,1) - bwb(k,:,2))*soil%zse ) +         &
+                          (bwb_gw(k,1) - bwb_gw(k,2))*soil%GWdz(k) ) *1000.0
+          END DO
+       ELSE IF(MOD(REAL(ktau),2.0)==0.0) THEN    ! if even timestep
+          bwb(:,:,2)=ssnow%wb
+          bwb_gw(:,2)=ssnow%GWwb
+          DO k=1,mp           !  current smoist - prev tstep smoist
+             delwb(k) = ( SUM( (bwb(k,:,2) - bwb(k,:,1))*soil%zse ) +         &
+                          (bwb_gw(k,2) -bwb_gw(k,1))*soil%GWdz(k) ) *1000.0
+          END DO
+       END IF
+      ! ______________________________________________________________________
+
+      ELSE ! MMY  ! IF part above inserted by rk4417 - phase2
+
        IF(MOD(REAL(ktau),2.0)==1.0) THEN         ! if odd timestep
           bwb(:,:,1)=ssnow%wb
           DO k=1,mp           ! current smoist - prev tstep smoist
@@ -513,7 +545,10 @@ SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met,
                   - (bwb(k,:,1)))*soil%zse)*1000.0
           END DO
        END IF
-    END IF
+
+    END IF ! MMY   ! inserted by rk4417 - phase2
+
+ END IF
 
 
 
@@ -523,10 +558,21 @@ SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met,
     !      it's included in change in canopy storage calculation))
     ! rml 28/2/11 ! BP changed rnof1+rnof2 to ssnow%runoff which also included rnof5
     ! which is used when nglacier=2 in soilsnow routines (BP feb2011)
+
+ IF ( cable_user%GW_MODEL) THEN ! MMY
+   ! ________________ MMY, Water Balance Equation for GW_ON _______________
+     bal%wbal = REAL(met%precip - canopy%delwc - ssnow%snowd+ssnow%osnowd       &
+         - ssnow%runoff-(canopy%fevw+canopy%fevc                                &
+         + canopy%fes/ssnow%cls)*dels/air%rlam - delwb) ! remove qrecharge
+   ! ______________________________________________________________________
+  ELSE ! MMY  ! IF part above inserted by rk4417 - phase2
+
     bal%wbal = REAL(met%precip - canopy%delwc - ssnow%snowd+ssnow%osnowd        &
          - ssnow%runoff-(canopy%fevw+canopy%fevc                                &
          + canopy%fes/ssnow%cls)*dels/air%rlam - delwb - ssnow%qrecharge)
 
+ END IF ! MMY   ! inserted by rk4417 - phase2
+
     ! Canopy water balance: precip-change.can.storage-throughfall-evap+dew
     canopy_wbal = REAL(met%precip-canopy%delwc-canopy%through                   &
          - (canopy%fevw+MIN(canopy%fevc,0.0_r_2))*dels/air%rlam)

src/offline/cable_define_types.F90

@@ -210,8 +210,9 @@ MODULE cable_def_types_mod
 
   ! Soil and snow variables:
   TYPE soil_snow_type
-
-     INTEGER, DIMENSION(:), POINTER :: isflag ! 0 => no snow 1 => snow
+     !> isflag 0 => one snow layer 1 => three snow layer
+     INTEGER, DIMENSION(:), POINTER :: isflag
+! 0 => one snow layer 1 => three snow layer ! MMY ! above and here inserted by rk4417 - phase2
 
      REAL, DIMENSION(:), POINTER ::                                           &
           iantrct, & ! pointer to Antarctic land points

src/offline/cable_driver.F90

@@ -62,7 +62,8 @@ PROGRAM cable_offline_driver
   USE cable_def_types_mod
   USE cable_IO_vars_module, ONLY: logn,gswpfile,ncciy,leaps,                  &
        verbose, fixedCO2,output,check,patchout,    &
-       patch_type,landpt,soilparmnew,&
+! line below modified by rk4417 - phase2
+       patch_type,landpt,& !soilparmnew,! MMY @Oct2022 change to use soilparmnew by default
        defaultLAI, sdoy, smoy, syear, timeunits, exists, calendar, set_group_output_values, &
        NO_CHECK
   USE casa_ncdf_module, ONLY: is_casa_time
@@ -245,7 +246,8 @@ PROGRAM cable_offline_driver
   NAMELIST/CABLE/                  &
        filename,         & ! TYPE, containing input filenames
        vegparmnew,       & ! use new soil param. method
-       soilparmnew,      & ! use new soil param. method
+!line below commented out by rk4417 - phase2
+!       soilparmnew,      & ! use new soil param. method ! MMY @Oct2022 change to use soilparmnew by default
        calcsoilalbedo,   & ! albedo considers soil color Ticket #27
        spinup,           & ! spinup model (soil) to steady state
        delsoilM,delsoilT,& !
@@ -409,8 +411,9 @@ PROGRAM cable_offline_driver
   !   STOP 'casaCNP required to get prognostic LAI or Vcmax'
   IF( l_vcmaxFeedbk .AND. icycle < 1 )                                     &
        STOP 'icycle must be 2 to 3 to get prognostic Vcmax'
-  IF( icycle > 0 .AND. ( .NOT. soilparmnew ) )                             &
-       STOP 'casaCNP must use new soil parameters'
+!line below commented out by rk4417 - phase2
+!  IF( icycle > 0 .AND. ( .NOT. soilparmnew ) ) & ! MMY @Oct2022 change to use soilparmnew by default
+!       STOP 'casaCNP must use new soil parameters'
 
   NRRRR = MERGE(MAX(CABLE_USER%CASA_NREP,1), 1, CASAONLY)
   ! casa time count
@@ -512,6 +515,40 @@ PROGRAM cable_offline_driver
                  kend      = ktauday * LOY
               ENDIF
 
+! block below inserted by rk4417 - phase2
+! __________________________ MMY using Princeton _______________________________
+        ELSE IF ( TRIM(cable_user%MetType) .EQ. 'prin' ) THEN
+           ncciy = CurYear
+
+           CALL prepareFiles_princeton(ncciy) ! MMY
+           IF ( RRRR .EQ. 1 ) THEN
+           CALL open_met_file( dels, koffset, kend, spinup, CTFRZ )
+           IF (leaps.and.is_leapyear(YYYY).and.kend.eq.2920) THEN
+              STOP 'LEAP YEAR INCOMPATIBILITY WITH INPUT MET !'
+           ENDIF
+           IF ( NRRRR .GT. 1 ) THEN
+              GSWP_MID(1,YYYY) = ncid_rain
+              ! GSWP_MID(2,YYYY) = ncid_snow MMY
+              GSWP_MID(3,YYYY) = ncid_lw
+              GSWP_MID(4,YYYY) = ncid_sw
+              GSWP_MID(5,YYYY) = ncid_ps
+              GSWP_MID(6,YYYY) = ncid_qa
+              GSWP_MID(7,YYYY) = ncid_ta
+              GSWP_MID(8,YYYY) = ncid_wd
+           ENDIF
+           ELSE
+           ncid_rain = GSWP_MID(1,YYYY)
+           ! ncid_snow = GSWP_MID(2,YYYY) MMY
+           ncid_lw   = GSWP_MID(3,YYYY)
+           ncid_sw   = GSWP_MID(4,YYYY)
+           ncid_ps   = GSWP_MID(5,YYYY)
+           ncid_qa   = GSWP_MID(6,YYYY)
+           ncid_ta   = GSWP_MID(7,YYYY)
+           ncid_wd   = GSWP_MID(8,YYYY)
+           kend      = ktauday * LOY ! MMY
+           ENDIF
+! ______________________________________________________________________________
+
            ELSE IF ( TRIM(cable_user%MetType) .EQ. 'plum' ) THEN
               ! PLUME experiment setup using WATCH
               IF ( CALL1 ) THEN
@@ -1049,9 +1086,17 @@ PROGRAM cable_offline_driver
                    YYYY.EQ. CABLE_USER%YearEnd ) THEN
 
                  ! evaluate spinup
+                 ! =================== MMY_phase2 commented out this region in favor of the one below - rk4417 =====================
                  IF( ANY( ABS(ssnow%wb-soilMtemp)>delsoilM).OR.               &
-                      ANY( ABS(ssnow%tgg-soilTtemp)>delsoilT) .OR. &
-                      MAXVAL(ABS(ssnow%GWwb-GWtemp),dim=1) > delgwM) THEN
+                 ANY( ABS(ssnow%tgg-soilTtemp)>delsoilT) .OR. &
+                 MAXVAL(ABS(ssnow%GWwb-GWtemp),dim=1) > delgwM) THEN
+
+                 ! =================== MMY_phase2 uncomment =====================                 
+!                 IF( (ANY( ABS(ssnow%wb-soilMtemp)>delsoilM).OR.                & 
+!                      ANY( ABS(ssnow%tgg-soilTtemp)>delsoilT) .or. &               
+!                      maxval(ABS(ssnow%GWwb-GWtemp),dim=1) > delgwM) .and. &
+!                      ( (int(ktau_tot/kend) .lt. cable_user%max_spins)  .and.&
+!                      (cable_user%max_spins .gt. 0) ) ) THEN
 
                     ! No complete convergence yet
                     PRINT *, 'ssnow%wb : ', ssnow%wb
@@ -1185,7 +1230,7 @@ PROGRAM cable_offline_driver
 
   IF (icycle > 0.and. .not.l_landuse) THEN
 
-     !CALL casa_poolout( ktau, veg, soil, casabiome,              &
+     !CALL casa_poolout( ktau, veg, soil, casabiome, & ! FEEDBACK (this call is uncommented in MMY code) --rk4417 
      ! casapool, casaflux, casamet, casabal, phen )
      CALL write_casa_restart_nc ( casamet, casapool,casaflux,phen, CASAONLY )
 
@@ -1307,6 +1352,50 @@ SUBROUTINE renameFiles(logn,inFile,ncciy,inName)
 
 END SUBROUTINE renameFiles
 
+
+! 2 subroutines below inserted by rk4417 - phase2
+! _______________________ MMY for Princeton input ______________________________
+SUBROUTINE prepareFiles_princeton(ncciy)
+  USE cable_IO_vars_module, ONLY: logn,gswpfile
+  IMPLICIT NONE
+  INTEGER, INTENT(IN) :: ncciy
+
+  WRITE(logn,*) 'CABLE offline global run using princeton forcing for ', ncciy
+  PRINT *,     'CABLE offline global run using princeton forcing for ', ncciy
+
+  CALL renameFiles_princeton(logn,gswpfile%rainf,ncciy,'rainf')
+  CALL renameFiles_princeton(logn,gswpfile%LWdown,ncciy,'LWdown')
+  CALL renameFiles_princeton(logn,gswpfile%SWdown,ncciy,'SWdown')
+  CALL renameFiles_princeton(logn,gswpfile%PSurf,ncciy,'PSurf')
+  CALL renameFiles_princeton(logn,gswpfile%Qair,ncciy,'Qair')
+  CALL renameFiles_princeton(logn,gswpfile%Tair,ncciy,'Tair')
+  CALL renameFiles_princeton(logn,gswpfile%wind,ncciy,'wind')
+
+END SUBROUTINE prepareFiles_princeton
+
+SUBROUTINE renameFiles_princeton(logn,inFile,ncciy,inName)
+  IMPLICIT NONE
+  INTEGER, INTENT(IN) :: logn,ncciy
+  INTEGER:: nn
+  CHARACTER(LEN=200), INTENT(INOUT) :: inFile
+  CHARACTER(LEN=*),  INTENT(IN)        :: inName
+  INTEGER :: idummy
+
+  nn = INDEX(inFile,'19')
+  READ(inFile(nn:nn+3),'(i4)') idummy
+  WRITE(inFile(nn:nn+3),'(i4.4)') ncciy
+  nn = INDEX(inFile,'19')
+  READ(inFile(nn:nn+3),'(i4)') idummy
+  WRITE(inFile(nn:nn+3),'(i4.4)') ncciy
+  READ(inFile(nn-5:nn-2),'(i4)') idummy
+  WRITE(inFile(nn-5:nn-2),'(i4.4)') ncciy
+  WRITE(logn,*) TRIM(inName), ' global data from ', TRIM(inFile)
+
+END SUBROUTINE renameFiles_princeton
+
+! ______________________________________________________________________________
+
+
 !==============================================================================
 ! subroutine for reading LU input data, zeroing biomass in empty secondary forest tiles
 ! and tranferring LUC-based age weights for secondary forest to POP structure