Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

477 add princeton met forcing #524

Draft
wants to merge 4 commits into
base: main
Choose a base branch
from
+216 −65
Draft

477 add princeton met forcing #524

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6213eaf
(#477) - Add princeton to input and driver
ccarouge Jan 13, 2025
5969248
No snow file for Princeton met forcing
ccarouge Jan 13, 2025
bbe15ff
Updates for new driver
ccarouge Jan 13, 2025
a960b32
(#525): Fix renameFiles_princeton
ccarouge Jan 14, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 48 additions & 3 deletions src/offline/cable_driver_common.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,7 @@ MODULE cable_driver_common_mod
PUBLIC :: cable_driver_init_default
PUBLIC :: prepareFiles
PUBLIC :: renameFiles
PUBLIC :: prepareFiles_princeton
PUBLIC :: LUCdriver

CONTAINS
Expand Down Expand Up @@ -332,7 +333,8 @@ SUBROUTINE cable_driver_init_cru(dels, koffset, CRU)
END SUBROUTINE cable_driver_init_cru

SUBROUTINE prepareFiles(ncciy)
!* Select the correct files given the year for filenames following the gswp format
!* Select the correct files given the year for filenames following
! the gswp format

USE cable_IO_vars_module, ONLY: logn,gswpfile
IMPLICIT NONE
Expand All @@ -350,10 +352,11 @@ SUBROUTINE prepareFiles(ncciy)
CALL renameFiles(logn,gswpfile%Tair,ncciy,'Tair')
CALL renameFiles(logn,gswpfile%wind,ncciy,'wind')

END SUBROUTINE prepareFiles
END SUBROUTINE prepareFiles

SUBROUTINE renameFiles(logn,inFile,ncciy,inName)
!! Replace the year in the filename with the value of ncciy.
!* Replace the year in the filename with the value of ncciy
! for the gswp file format.

IMPLICIT NONE
INTEGER, INTENT(IN) :: logn !! Log file unit number
Expand All @@ -370,6 +373,48 @@ SUBROUTINE renameFiles(logn,inFile,ncciy,inName)

END SUBROUTINE renameFiles

SUBROUTINE prepareFiles_princeton(ncciy)
!* Select the correct files given the year for filenames following the
! princeton format
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)
!* Replace the year in the filename with the value of ncciy for
! the princeton format
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', BACK=.TRUE.)
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
SUBROUTINE LUCdriver( casabiome,casapool, &
Expand Down
207 changes: 149 additions & 58 deletions src/offline/cable_input.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -400,7 +400,7 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ)
PRINT*,'rainf'
CALL handle_err( ok )
ENDIF
IF(.NOT. globalMetfile%l_gpcc) THEN
IF(.NOT. globalMetfile%l_gpcc .AND. cable_user%MetType .NE. "prin") THEN
ok = NF90_OPEN(gswpfile%snowf,0,ncid_snow)
IF (ok /= NF90_NOERR) THEN
PRINT*,'snow'
Expand Down 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
11 changes: 9 additions & 2 deletions src/offline/cable_mpimaster.F90
View file
Open in desktop
Copy link
Member Author

@ccarouge ccarouge Jan 16, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ccarouge might need something more. See l.877 in the GW PR.
There should be something for the output as well, see l. 989 in the GW PR.

Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,7 @@ MODULE cable_mpimaster
LALLOC, &
prepareFiles, &
renameFiles, &
prepareFiles_princeton, &
LUCdriver
USE cable_mpicommon
USE cable_IO_vars_module, ONLY : NO_CHECK
Expand Down Expand Up @@ -353,11 +354,17 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
ENDIF

SELECT CASE (TRIM(cable_user%MetType))
CASE ('gswp')
CASE ('gswp', 'prin')
ncciy = CurYear

WRITE(*,*) 'Looking for global offline run info.'
CALL prepareFiles(ncciy)

IF ( TRIM(cable_user%MetType) == 'gswp' ) THEN
CALL prepareFiles(ncciy)
ELSE ! cable_user%MetType == 'prin'
CALL prepareFiles_princeton(ncciy)
END IF

CALL open_met_file( dels, koffset, kend, spinup, CTFRZ )

CASE ('plum')
Expand Down
Loading
Loading