Skip to content

Commit

Permalink
(#477) - Add princeton to input and driver
Browse files Browse the repository at this point in the history
  • Loading branch information
ccarouge committed Dec 3, 2024
1 parent dfef1b9 commit de1be55
Show file tree
Hide file tree
Showing 2 changed files with 225 additions and 58 deletions.
205 changes: 148 additions & 57 deletions src/offline/cable_input.F90
Original file line number Diff line number Diff line change
Expand Up @@ -497,9 +497,12 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ)
IF(ok/=NF90_NOERR) THEN ! if failed
! Try 'lon' instead of x
ok = NF90_INQ_DIMID(ncid_met,'lon', xdimID)
IF(ok/=NF90_NOERR) CALL nc_abort &
(ok,'Error finding x dimension in '&
//TRIM(filename%met)//' (SUBROUTINE open_met_file)')
IF(ok/=NF90_NOERR) THEN ! MMY
ok = NF90_INQ_DIMID(ncid_met,'longitude', xdimID) ! MMY ! For princeton
IF(ok/=NF90_NOERR) CALL nc_abort & ! MMY
(ok,'Error finding x dimension in '& ! MMY
//TRIM(filename%met)//' (SUBROUTINE open_met_file)') ! MMY
END IF ! MMY
END IF
ok = NF90_INQUIRE_DIMENSION(ncid_met,xdimID,len=xdimsize)
IF(ok/=NF90_NOERR) CALL nc_abort &
Expand All @@ -510,9 +513,12 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ)
IF(ok/=NF90_NOERR) THEN ! if failed
! Try 'lat' instead of y
ok = NF90_INQ_DIMID(ncid_met,'lat', ydimID)
IF(ok/=NF90_NOERR) CALL nc_abort &
(ok,'Error finding y dimension in ' &
//TRIM(filename%met)//' (SUBROUTINE open_met_file)')
IF(ok/=NF90_NOERR) THEN ! MMY
ok = NF90_INQ_DIMID(ncid_met,'latitude', ydimID) ! MMY ! For princeton
IF(ok/=NF90_NOERR) CALL nc_abort & ! MMY
(ok,'Error finding y dimension in ' & ! MMY
//TRIM(filename%met)//' (SUBROUTINE open_met_file)') ! MMY
END IF ! MMY
END IF
ok = NF90_INQUIRE_DIMENSION(ncid_met,ydimID,len=ydimsize)
IF(ok/=NF90_NOERR) CALL nc_abort &
Expand Down Expand Up @@ -821,7 +827,10 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ)
!===== done bug fixing for timevar in PALS met file ===============

!===== gswp input file has bug in timeunits ===========
IF (ncciy > 0) WRITE(timeunits(26:27),'(i2.2)') 0
IF (TRIM(cable_user%MetType) .NE. "prin") THEN ! MMY
IF (ncciy > 0) WRITE(timeunits(26:27),'(i2.2)') 0
END IF ! MMY

!===== done bug fixing for timeunits in gwsp file ========
WRITE(logn,*) 'Time variable units: ', timeunits
! Get coordinate field:
Expand All @@ -848,17 +857,25 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ)
! start time) from character to integer; calculate starting hour-of-day,
! day-of-year, year:
IF (.NOT.cable_user%GSWP3) THEN
READ(timeunits(15:18),*) syear
READ(timeunits(20:21),*) smoy ! integer month
READ(timeunits(23:24),*) sdoytmp ! integer day of that month
READ(timeunits(26:27),*) shod ! starting hour of day
IF (cable_user%MetType .eq. "prin") THEN ! MMY
READ(timeunits(13:16),*) syear ! MMY
READ(timeunits(18:19),*) smoy ! integer month ! MMY
READ(timeunits(21:22),*) sdoytmp ! integer day of that month ! MMY
READ(timeunits(24:25),*) shod ! starting hour of day ! MMY
ELSE ! MMY
READ(timeunits(15:18),*) syear
READ(timeunits(20:21),*) smoy ! integer month
READ(timeunits(23:24),*) sdoytmp ! integer day of that month
READ(timeunits(26:27),*) shod ! starting hour of day
END IF
ELSE
syear=ncciy
smoy=1
sdoytmp=1
shod=0
END IF


! if site data, shift start time to middle of timestep
! only do this if not already at middle of timestep
! vh_js !
Expand Down Expand Up @@ -1055,14 +1072,19 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ)
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error finding Rainf units in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE open_met_file)')
IF(metunits%Rainf(1:8)=='kg/m^2/s'.OR.metunits%Rainf(1:6)=='kg/m2s'.OR.metunits%Rainf(1:10)== &
'kgm^-2s^-1'.OR.metunits%Rainf(1:4)=='mm/s'.OR. &
metunits%Rainf(1:6)=='mms^-1'.OR. &
metunits%Rainf(1:7)=='kg/m^2s'.OR.metunits%Rainf(1:10)=='kg m-2 s-1'.OR.metunits%Wind(1:5)/='m s-1') THEN
IF(metunits%Rainf(1:8)=='kg/m^2/s' .OR. &
metunits%Rainf(1:7)=='kg/m2/s' .OR. &
metunits%Rainf(1:6)=='kg/m2s' .OR. &
metunits%Rainf(1:10)=='kgm^-2s^-1' .OR. &
metunits%Rainf(1:4)=='mm/s' .OR. &
metunits%Rainf(1:6)=='mms^-1' .OR. &
metunits%Rainf(1:7)=='kg/m^2s' .OR. &
metunits%Rainf(1:10)=='kg m-2 s-1' .OR. &
metunits%Wind(1:5)/='m s-1' ) THEN
! Change from mm/s to mm/time step:
convert%Rainf = dels
ELSE IF(metunits%Rainf(1:4)=='mm/h'.OR.metunits%Rainf(1:6)== &
'mmh^-1') THEN
ELSE IF(metunits%Rainf(1:4)=='mm/h' .OR. &
metunits%Rainf(1:6)== 'mmh^-1' ) THEN
! Change from mm/h to mm/time step:
convert%Rainf = dels/3600.0
ELSE
Expand Down Expand Up @@ -1817,19 +1839,39 @@ SUBROUTINE get_met_data(spinup,spinConv,met,soil,rad, &

! Get SWdown data for mask grid:
IF (cable_user%GSWP3) ncid_met=ncid_sw ! since GSWP3 multiple met files
ok= NF90_GET_VAR(ncid_met,id%SWdown,tmpDat3, &
start=(/1,1,ktau/),count=(/xdimsize,ydimsize,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading SWdown in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
! Assign value to met data variable (no units change required):
!jhan:quick fix, use (1/2) dimension 1 here arbitrarily
DO i=1,mland ! over all land points/grid cells
met%fsd(landpt(i)%cstart:landpt(i)%cend,1) = &
0.5 * REAL(tmpDat3(land_x(i),land_y(i),1))
met%fsd(landpt(i)%cstart:landpt(i)%cend,2) = &
0.5 * REAL(tmpDat3(land_x(i),land_y(i),1))
ENDDO

! ______________________ MMY _________________________
IF (TRIM(cable_user%MetType) .EQ. "prin") THEN
ok= NF90_GET_VAR(ncid_met,id%SWdown,tmpDat4, &
start=(/1,1,1,ktau/),count=(/xdimsize,ydimsize,1,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading SWdown in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%fsd(landpt(i)%cstart:landpt(i)%cend,1) = &
0.5 * REAL(tmpDat4(land_x(i),land_y(i),1,1))
met%fsd(landpt(i)%cstart:landpt(i)%cend,2) = &
0.5 * REAL(tmpDat4(land_x(i),land_y(i),1,1))
ENDDO
! PRINT *, "========== MMY =========="
! PRINT *, "met%fsd",met%fsd
! ____________________________________________________
ELSE ! MMY

ok= NF90_GET_VAR(ncid_met,id%SWdown,tmpDat3, &
start=(/1,1,ktau/),count=(/xdimsize,ydimsize,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading SWdown in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
! Assign value to met data variable (no units change required):
!jhan:quick fix, use (1/2) dimension 1 here arbitrarily
DO i=1,mland ! over all land points/grid cells
met%fsd(landpt(i)%cstart:landpt(i)%cend,1) = &
0.5 * REAL(tmpDat3(land_x(i),land_y(i),1))
met%fsd(landpt(i)%cstart:landpt(i)%cend,2) = &
0.5 * REAL(tmpDat3(land_x(i),land_y(i),1))
ENDDO
END IF ! MMY ! inserted by rk4417 - phase2

! Get Tair data for mask grid:- - - - - - - - - - - - - - - - - -
IF(cable_user%GSWP3) ncid_met = ncid_ta ! since GSWP3 multiple met files
Expand Down Expand Up @@ -2022,15 +2064,31 @@ SUBROUTINE get_met_data(spinup,spinConv,met,soil,rad, &

! Get Rainf and Snowf data for mask grid:- - - - - - - - - - - - -
IF (cable_user%GSWP3) ncid_met = ncid_rain
ok= NF90_GET_VAR(ncid_met,id%Rainf,tmpDat3, &
start=(/1,1,ktau/),count=(/xdimsize,ydimsize,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading Rainf in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%precip(landpt(i)%cstart:landpt(i)%cend) = &
REAL(tmpDat3(land_x(i),land_y(i),1)) ! store Rainf
ENDDO

! ______________________ MMY _________________________
IF (TRIM(cable_user%MetType) .eq. "prin") THEN ! MMY
ok= NF90_GET_VAR(ncid_met,id%Rainf,tmpDat4, &
start=(/1,1,1,ktau/),count=(/xdimsize,ydimsize,1,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading Rainf in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%precip(landpt(i)%cstart:landpt(i)%cend) = &
REAL(tmpDat4(land_x(i),land_y(i),1,1)) ! store Rainf
ENDDO
! ____________________________________________________
ELSE ! MMY
ok= NF90_GET_VAR(ncid_met,id%Rainf,tmpDat3, &
start=(/1,1,ktau/),count=(/xdimsize,ydimsize,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading Rainf in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%precip(landpt(i)%cstart:landpt(i)%cend) = &
REAL(tmpDat3(land_x(i),land_y(i),1)) ! store Rainf
ENDDO
END IF

IF(exists%Snowf) THEN
IF (cable_user%GSWP3) ncid_met = ncid_snow
ok= NF90_GET_VAR(ncid_met,id%Snowf,tmpDat3, &
Expand Down Expand Up @@ -2067,32 +2125,65 @@ SUBROUTINE get_met_data(spinup,spinConv,met,soil,rad, &

! Get LWdown data for mask grid: - - - - - - - - - - - - - - - - -
IF (cable_user%GSWP3) ncid_met = ncid_lw
IF(exists%LWdown) THEN ! If LWdown exists in met file
ok= NF90_GET_VAR(ncid_met,id%LWdown,tmpDat3, &
start=(/1,1,ktau/),count=(/xdimsize,ydimsize,1/))

! ______________________ MMY _________________________
IF (TRIM(cable_user%MetType) .eq. "prin") THEN
ok= NF90_GET_VAR(ncid_met,id%LWdown,tmpDat4, &
start=(/1,1,1,ktau/),count=(/xdimsize,ydimsize,1,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading LWdown in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
(ok,'Error reading LWdown in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%fld(landpt(i)%cstart:landpt(i)%cend) = &
REAL(tmpDat3(land_x(i),land_y(i),1))
REAL(tmpDat4(land_x(i),land_y(i),1,1))
ENDDO
ELSE ! Synthesise LWdown based on temperature
! Use Swinbank formula:
met%fld(:) = 0.0000094*0.0000000567*(met%tk(:)**6.0)
! PRINT *, "========== MMY =========="
! PRINT *, "met%fld",met%fld
! ____________________________________________________
ELSE ! MMY

IF(exists%LWdown) THEN ! If LWdown exists in met file
ok= NF90_GET_VAR(ncid_met,id%LWdown,tmpDat3, &
start=(/1,1,ktau/),count=(/xdimsize,ydimsize,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading LWdown in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%fld(landpt(i)%cstart:landpt(i)%cend) = &
REAL(tmpDat3(land_x(i),land_y(i),1))
ENDDO
ELSE ! Synthesise LWdown based on temperature
! Use Swinbank formula:
met%fld(:) = 0.0000094*0.0000000567*(met%tk(:)**6.0)
END IF
END IF

! Get CO2air data for mask grid:- - - - - - - - - - - - - - - - - -
IF(exists%CO2air) THEN ! If CO2air exists in met file
ok= NF90_GET_VAR(ncid_met,id%CO2air,tmpDat4, &
start=(/1,1,1,ktau/),count=(/xdimsize,ydimsize,1,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading CO2air in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%ca(landpt(i)%cstart:landpt(i)%cend) = &
REAL(tmpDat4(land_x(i),land_y(i),1,1))/1000000.0
ENDDO
! ____________ MMY@23Apr2023 to read CO2air from PLUMBER2 __________
ok = NF90_INQUIRE_VARIABLE(ncid_met,id%CO2air,ndims=ndims)
IF(ndims==3) THEN
ok= NF90_GET_VAR(ncid_met,id%CO2air,tmpDat3, &
start=(/1,1,ktau/),count=(/xdimsize,ydimsize,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading CO2air in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%ca(landpt(i)%cstart:landpt(i)%cend) = &
REAL(tmpDat3(land_x(i),land_y(i),1))/1000000.0
ENDDO
ELSE
! __________________________________________________________________
ok= NF90_GET_VAR(ncid_met,id%CO2air,tmpDat4, &
start=(/1,1,1,ktau/),count=(/xdimsize,ydimsize,1,1/))
IF(ok /= NF90_NOERR) CALL nc_abort &
(ok,'Error reading CO2air in met data file ' &
//TRIM(filename%met)//' (SUBROUTINE get_met_data)')
DO i=1,mland ! over all land points/grid cells
met%ca(landpt(i)%cstart:landpt(i)%cend) = &
REAL(tmpDat4(land_x(i),land_y(i),1,1))/1000000.0
ENDDO
END IF
ELSE
! Fix CO2 air concentration:
met%ca(:) = fixedCO2 /1000000.0
Expand Down
78 changes: 77 additions & 1 deletion src/offline/cable_serial.F90
Original file line number Diff line number Diff line change
Expand Up @@ -322,7 +322,39 @@ SUBROUTINE serialdrv(trunk_sumbal, NRRRR, dels, koffset, kend, GSWP_MID, PLUME,
kend = ktauday * LOY
ENDIF

ELSE IF ( TRIM(cable_user%MetType) .EQ. 'plum' ) THEN
! __________________________ 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 ( .NOT. PLUME%LeapYears ) LOY = 365
kend = NINT(24.0*3600.0/dels) * LOY
Expand Down Expand Up @@ -1063,6 +1095,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
Expand Down

0 comments on commit de1be55

Please sign in to comment.