From de1be55fe1712f68bc52a144bbdcf60a4c3de36a Mon Sep 17 00:00:00 2001 From: Claire Carouge Date: Tue, 3 Dec 2024 11:01:28 +1100 Subject: [PATCH] (#477) - Add princeton to input and driver --- src/offline/cable_input.F90 | 205 +++++++++++++++++++++++++---------- src/offline/cable_serial.F90 | 78 ++++++++++++- 2 files changed, 225 insertions(+), 58 deletions(-) diff --git a/src/offline/cable_input.F90 b/src/offline/cable_input.F90 index 9c30208f7..84d444db7 100644 --- a/src/offline/cable_input.F90 +++ b/src/offline/cable_input.F90 @@ -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 & @@ -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 & @@ -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: @@ -848,10 +857,17 @@ 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 @@ -859,6 +875,7 @@ SUBROUTINE open_met_file(dels,koffset,kend,spinup, TFRZ) 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 ! @@ -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 @@ -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 @@ -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, & @@ -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 diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90 index 11df02a73..c3f66ddbc 100644 --- a/src/offline/cable_serial.F90 +++ b/src/offline/cable_serial.F90 @@ -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 @@ -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