This code is provided as is by Dag Lohmann
The include files needed by this program are :
Instruction to compile the code on the NCEP system are provided in the Makefile
IMPLICIT NONE
include '/ldas/home/noah/PILPS-2E/netcdf-3.4/include/netcdf.inc'
INCLUDE 'MAIN.h'
include 'netcdf2.h'
C DRIVER STEP 1
C READ CONTROL FILE
WRITE(*,*) 'READ CONTROLFILE'
CNTRFL = 'controlfile'
CALL READCNTL(CNTRFL,ICE,DT,Z,SNOALB)
C DRIVER STEP 2
C READ MODEL GEOMETRY / LAND SURFACE CHARACTERISTICS
WRITE(*,*) 'READ MODEL GEOMETRY'
CALL READ_LSC(NX,NY,NSOLD,NMONTH,NSOIL,SOILTYP,
& LAND_SEA,VEGTYP,SLOPETYP,SOILDEPTH,ALBEDO,SHDFAC,
& TBOT)
C DRIVER STEP 3
C READ INITIAL CONDITIONS
WRITE(*,*) 'READ INITIAL CONDITIONS'
CALL READ_INITIAL(NX,NY,NSOLD,T1,STC,SMC,SH2O,CMC,
& SNOWH,SNEQV,CH,CM)
C CHECK INITIAL CONDITIONS AND MODEL GEOMETRY
CALL CHECK_INITIAL(NX,NY,NSOLD,LAND_SEA,T1,STC,SMC,SH2O,CMC,
& SNOWH,SNEQV,SOILTYP,VEGTYP,SLOPETYP,TBOT,ALBEDO,SHDFAC,
& NSOIL,NMONTH,SOILDEPTH,CH,CM)
C DRIVER STEP 5
C RUN TIME AND SPATIAL LOOP
C timestp = 87672
timestp = 0
DO IYEAR = 1,20
FIRST_STEP = .TRUE.
WRITE(YEAR_CH,'(I2.2)') IYEAR+9
WRITE(YEAR_CH2,'(I2.2)') IYEAR+78
FILENAME = 'fort.'//YEAR_CH//char(0)
WRITE(*,*) 'reading atmospheric data ', FILENAME
flag = 0
status = NF_OPEN(FILENAME, flag, ncid)
status = NF_INQ(ncid, ndims, nvars, ngatts, unlimited)
status = NF_INQ_DIMID(ncid,'x',xid)
status = NF_INQ_DIMID(ncid,'y',yid)
status = NF_INQ_DIMID(ncid,'tstep',tstepid)
status = NF_INQ_DIMLEN(ncid,xid,xlen)
status = NF_INQ_DIMLEN(ncid,yid,ylen)
status = NF_INQ_DIMLEN(ncid,tstepid,tsteplen)
C CREATE OUTPUT FILE FOR EACH YEAR
out1name = 'noah.rerun.eb.torne.hrly.19'//YEAR_CH2//'.nc'
flag = 0
status = NF_CREATE(out1name,flag,out1)
status = NF_DEF_DIM(out1,'x',xlen,xid)
status = NF_DEF_DIM(out1,'y',ylen,yid)
status = NF_DEF_DIM(out1,'tstep',NF_unlimited,tstepid)
outdim(1) = xid
outdim(2) = yid
outdim(3) = tstepid
status = NF_DEF_VAR(out1,'lon',NF_FLOAT,1,xid,lonid1)
status = NF_DEF_VAR(out1,'lat',NF_FLOAT,1,yid,latid1)
status = NF_DEF_VAR(out1,'time',NF_FLOAT,1,tstepid,timeid1)
status = NF_DEF_VAR(out1,'timestp',NF_INT,1,tstepid,timestpid1)
status = NF_PUT_ATT_TEXT(out1,lonid1,'units',12,'degrees_east')
status = NF_PUT_ATT_REAL(out1,lonid1,'valid_min',NF_FLOAT,1,
& 18.125)
status = NF_PUT_ATT_REAL(out1,lonid1,'valid_max',NF_FLOAT,1,
& 25.125)
status = NF_PUT_ATT_TEXT(out1,lonid1,'long_name',9,'Longitude')
status = NF_PUT_ATT_TEXT(out1,latid1,'units',13,
& 'degrees_north')
status = NF_PUT_ATT_REAL(out1,latid1,'valid_min',NF_FLOAT,1,
& 65.875)
status = NF_PUT_ATT_REAL(out1,latid1,'valid_max',NF_FLOAT,1,
& 69.125)
status = NF_PUT_ATT_TEXT(out1,latid1,'long_name',8,'Latitude')
status = NF_PUT_ATT_TEXT(out1,timeid1,'units',33,
& 'seconds since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out1,timeid1,'calendar',9,
& 'gregorian')
status = NF_PUT_ATT_TEXT(out1,timeid1,'title',4,
& 'Time')
status = NF_PUT_ATT_TEXT(out1,timeid1,'long_name',9,
& 'Time axis')
status = NF_PUT_ATT_TEXT(out1,timeid1,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_PUT_ATT_TEXT(out1,timestpid1,'units',35,
& 'timesteps since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out1,timestpid1,'title',10,
& 'Time steps')
status = NF_PUT_ATT_REAL(out1,timestpid1,'tstep_sec',
& NF_FLOAT,1,3600.0)
status = NF_PUT_ATT_TEXT(out1,timestpid1,'long_name',14,
& 'Time step axis')
status = NF_PUT_ATT_TEXT(out1,timestpid1,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_DEF_VAR(out1,'SWnet',NF_FLOAT,3,outdim,
& SWnetid)
status = NF_DEF_VAR(out1,'LWnet',NF_FLOAT,3,outdim,
& LWnetid)
status = NF_DEF_VAR(out1,'Qle',NF_FLOAT,3,outdim,Qleid)
status = NF_DEF_VAR(out1,'Qh',NF_FLOAT,3,outdim,Qhid)
status = NF_DEF_VAR(out1,'Qg',NF_FLOAT,3,outdim,Qgid)
status = NF_PUT_ATT_TEXT(out1,NF_GLOBAL,'Conventions',7,
& 'GDT 1.2')
status = NF_PUT_ATT_TEXT(out1,NF_GLOBAL,'file_name',32,
& 'noah.rerun.eb.torne.hrly.19'//YEAR_CH2//'.nc')
status = NF_PUT_ATT_TEXT(out1,NF_GLOBAL,'production',30,
& 'NOAH model run, DL on 20010312')
status = NF_PUT_ATT_TEXT(out1,SWnetid,'units',5,'W/m^2')
status = NF_PUT_ATT_TEXT(out1,SWnetid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out1,SWnetid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out1,SWnetid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out1,LWnetid,'units',5,'W/m^2')
status = NF_PUT_ATT_TEXT(out1,LWnetid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out1,LWnetid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out1,LWnetid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out1,Qleid,'units',5,'W/m^2')
status = NF_PUT_ATT_TEXT(out1,Qleid,'axis',11,'Fortran XYT')
status = NF_PUT_ATT_REAL(out1,Qleid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out1,Qleid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out1,Qhid,'units',5,'W/m^2')
status = NF_PUT_ATT_TEXT(out1,Qhid,'axis',11,'Fortran XYT')
status = NF_PUT_ATT_REAL(out1,Qhid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out1,Qhid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out1,Qgid,'units',5,'W/m^2')
status = NF_PUT_ATT_TEXT(out1,Qgid,'axis',11,'Fortran XYT')
status = NF_PUT_ATT_REAL(out1,Qgid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out1,Qgid,'associate',12,
& 'lon lat time')
status = NF_ENDDEF(out1)
out2name = 'noah.rerun.wb.torne.hrly.19'//YEAR_CH2//'.nc'
flag = 0
status = NF_CREATE(out2name,flag,out2)
status = NF_DEF_DIM(out2,'x',xlen,xid)
status = NF_DEF_DIM(out2,'y',ylen,yid)
status = NF_DEF_DIM(out2,'tstep',NF_unlimited,tstepid)
outdim(1) = xid
outdim(2) = yid
outdim(3) = tstepid
status = NF_DEF_VAR(out2,'lon',NF_FLOAT,1,xid,lonid2)
status = NF_DEF_VAR(out2,'lat',NF_FLOAT,1,yid,latid2)
status = NF_DEF_VAR(out2,'time',NF_FLOAT,1,tstepid,timeid2)
status = NF_DEF_VAR(out2,'timestp',NF_INT,1,tstepid,timestpid2)
status = NF_PUT_ATT_TEXT(out2,lonid2,'units',12,'degrees_east')
status = NF_PUT_ATT_REAL(out2,lonid2,'valid_min',NF_FLOAT,1,
& 18.125)
status = NF_PUT_ATT_REAL(out2,lonid2,'valid_max',NF_FLOAT,1,
& 25.125)
status = NF_PUT_ATT_TEXT(out2,lonid2,'long_name',9,'Longitude')
status = NF_PUT_ATT_TEXT(out2,latid2,'units',13,
& 'degrees_north')
status = NF_PUT_ATT_REAL(out2,latid2,'valid_min',NF_FLOAT,1,
& 65.875)
status = NF_PUT_ATT_REAL(out2,latid2,'valid_max',NF_FLOAT,1,
& 69.125)
status = NF_PUT_ATT_TEXT(out2,latid2,'long_name',8,'Latitude')
status = NF_PUT_ATT_TEXT(out2,timeid2,'units',33,
& 'seconds since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out2,timeid2,'calendar',9,
& 'gregorian')
status = NF_PUT_ATT_TEXT(out2,timeid2,'title',4,
& 'Time')
status = NF_PUT_ATT_TEXT(out2,timeid2,'long_name',9,
& 'Time axis')
status = NF_PUT_ATT_TEXT(out2,timeid2,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_PUT_ATT_TEXT(out2,timestpid2,'units',35,
& 'timesteps since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out2,timestpid2,'title',10,
& 'Time steps')
status = NF_PUT_ATT_REAL(out2,timestpid2,'tstep_sec',
& NF_FLOAT,1,3600.0)
status = NF_PUT_ATT_TEXT(out2,timestpid2,'long_name',14,
& 'Time step axis')
status = NF_PUT_ATT_TEXT(out2,timestpid2,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_DEF_VAR(out2,'Snowf',NF_FLOAT,3,outdim,
& Snowf2id)
status = NF_DEF_VAR(out2,'Rainf',NF_FLOAT,3,outdim,
& Rainf2id)
status = NF_DEF_VAR(out2,'Evap',NF_FLOAT,3,outdim,Evapid)
status = NF_DEF_VAR(out2,'Qs',NF_FLOAT,3,outdim,Qsid)
status = NF_DEF_VAR(out2,'Qsb',NF_FLOAT,3,outdim,Qsbid)
status = NF_DEF_VAR(out2,'Qsm',NF_FLOAT,3,outdim,Qsmid)
status = NF_DEF_VAR(out2,'DelSoilMoist',NF_FLOAT,3,outdim,
& DelSoilid)
status = NF_DEF_VAR(out2,'DelSWE',NF_FLOAT,3,outdim,DelSWEid)
status = NF_DEF_VAR(out2,'DelIntercept',NF_FLOAT,3,outdim,
& DelInterid)
status = NF_PUT_ATT_TEXT(out2,NF_GLOBAL,'Conventions',7,
& 'GDT 1.2')
status = NF_PUT_ATT_TEXT(out2,NF_GLOBAL,'file_name',32,
& 'noah.rerun.wb.torne.hrly.19'//YEAR_CH2//'.nc')
status = NF_PUT_ATT_TEXT(out2,NF_GLOBAL,'production',30,
& 'NOAH model run, DL on 20010312')
status = NF_PUT_ATT_TEXT(out2,Snowf2id,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out2,Snowf2id,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,Snowf2id,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,Snowf2id,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out2,Rainf2id,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out2,Rainf2id,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,Rainf2id,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,Rainf2id,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out2,Evapid,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out2,Evapid,'axis',11,'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,Evapid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,Evapid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out2,Qsid,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out2,Qsid,'axis',11,'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,Qsid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,Qsid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out2,Qsbid,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out2,Qsbid,'axis',11,'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,Qsbid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,Qsbid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out2,Qsmid,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out2,Qsmid,'axis',11,'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,Qsmid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,Qsmid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out2,DelSoilid,'units',6,'kg/m^2')
status = NF_PUT_ATT_TEXT(out2,DelSoilid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,DelSoilid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,DelSoilid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out2,DelSWEid,'units',6,'kg/m^2')
status = NF_PUT_ATT_TEXT(out2,DelSWEid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,DelSWEid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,DelSWEid,'associate',12,
& 'lon lat time')
status = NF_PUT_ATT_TEXT(out2,DelInterid,'units',6,'kg/m^2')
status = NF_PUT_ATT_TEXT(out2,DelInterid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out2,DelInterid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out2,DelInterid,'associate',12,
& 'lon lat time')
status = NF_ENDDEF(out2)
out3name = 'noah.rerun.sur.torne.hrly.19'//YEAR_CH2//'.nc'
flag = 0
status = NF_CREATE(out3name,flag,out3)
status = NF_DEF_DIM(out3,'x',xlen,xid)
status = NF_DEF_DIM(out3,'y',ylen,yid)
status = NF_DEF_DIM(out3,'tstep',NF_unlimited,tstepid)
outdim(1) = xid
outdim(2) = yid
outdim(3) = tstepid
status = NF_PUT_ATT_TEXT(out3,NF_GLOBAL,'Conventions',7,
& 'GDT 1.2')
status = NF_PUT_ATT_TEXT(out3,NF_GLOBAL,'file_name',33,
& 'noah.rerun.sur.torne.hrly.19'//YEAR_CH2//'.nc')
status = NF_PUT_ATT_TEXT(out3,NF_GLOBAL,'production',30,
& 'NOAH model run, DL on 20010312')
status = NF_DEF_VAR(out3,'lon',NF_FLOAT,1,xid,lonid3)
status = NF_DEF_VAR(out3,'lat',NF_FLOAT,1,yid,latid3)
status = NF_DEF_VAR(out3,'time',NF_FLOAT,1,tstepid,timeid3)
status = NF_DEF_VAR(out3,'timestp',NF_INT,1,tstepid,timestpid3)
status = NF_PUT_ATT_TEXT(out3,lonid3,'units',12,'degrees_east')
status = NF_PUT_ATT_REAL(out3,lonid3,'valid_min',NF_FLOAT,1,
& 18.125)
status = NF_PUT_ATT_REAL(out3,lonid3,'valid_max',NF_FLOAT,1,
& 25.125)
status = NF_PUT_ATT_TEXT(out3,lonid3,'long_name',9,'Longitude')
status = NF_PUT_ATT_TEXT(out3,latid3,'units',13,
& 'degrees_north')
status = NF_PUT_ATT_REAL(out3,latid3,'valid_min',NF_FLOAT,1,
& 65.875)
status = NF_PUT_ATT_REAL(out3,latid3,'valid_max',NF_FLOAT,1,
& 69.125)
status = NF_PUT_ATT_TEXT(out3,latid3,'long_name',8,'Latitude')
status = NF_PUT_ATT_TEXT(out3,timeid3,'units',33,
& 'seconds since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out3,timeid3,'calendar',9,
& 'gregorian')
status = NF_PUT_ATT_TEXT(out3,timeid3,'title',4,
& 'Time')
status = NF_PUT_ATT_TEXT(out3,timeid3,'long_name',9,
& 'Time axis')
status = NF_PUT_ATT_TEXT(out3,timeid3,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_PUT_ATT_TEXT(out3,timestpid3,'units',35,
& 'timesteps since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out3,timestpid3,'title',10,
& 'Time steps')
status = NF_PUT_ATT_REAL(out3,timestpid3,'tstep_sec',
& NF_FLOAT,1,3600.0)
status = NF_PUT_ATT_TEXT(out3,timestpid3,'long_name',14,
& 'Time step axis')
status = NF_PUT_ATT_TEXT(out3,timestpid3,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_DEF_VAR(out3,'SnowT',NF_FLOAT,3,outdim,
& SnowTid)
status = NF_PUT_ATT_TEXT(out3,SnowTid,'units',1,'K')
status = NF_PUT_ATT_TEXT(out3,SnowTid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out3,SnowTid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out3,SnowTid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out3,'BaresoilT',NF_FLOAT,3,outdim,
& BaresoilTid)
status = NF_PUT_ATT_TEXT(out3,BaresoilTid,'units',1,'K')
status = NF_PUT_ATT_TEXT(out3,BaresoilTid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out3,BaresoilTid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out3,BaresoilTid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out3,'RadT',NF_FLOAT,3,outdim,
& RadTid)
status = NF_PUT_ATT_TEXT(out3,RadTid,'units',1,'K')
status = NF_PUT_ATT_TEXT(out3,RadTid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out3,RadTid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out3,RadTid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out3,'SWE',NF_FLOAT,3,outdim,
& SWEid)
status = NF_PUT_ATT_TEXT(out3,SWEid,'units',6,'kg/m^2')
status = NF_PUT_ATT_TEXT(out3,SWEid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out3,SWEid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out3,SWEid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out3,'VegT',NF_FLOAT,3,outdim,
& VegTid)
status = NF_PUT_ATT_TEXT(out3,VegTid,'units',1,'K')
status = NF_PUT_ATT_TEXT(out3,VegTid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out3,VegTid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out3,VegTid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out3,'AvgSurfT',NF_FLOAT,3,outdim,
& AvgSurfTid)
status = NF_PUT_ATT_TEXT(out3,AvgSurfTid,'units',1,'K')
status = NF_PUT_ATT_TEXT(out3,AvgSurfTid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out3,AvgSurfTid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out3,AvgSurfTid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out3,'Albedo',NF_FLOAT,3,outdim,
& Albedoid)
status = NF_PUT_ATT_TEXT(out3,Albedoid,'units',1,'-')
status = NF_PUT_ATT_TEXT(out3,Albedoid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out3,Albedoid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out3,Albedoid,'associate',12,
& 'lon lat time')
status = NF_ENDDEF(out3)
out4name = 'noah.rerun.sub.torne.hrly.19'//YEAR_CH2//'.nc'
flag = 0
status = NF_CREATE(out4name,flag,out4)
status = NF_DEF_DIM(out4,'x',xlen,xid)
status = NF_DEF_DIM(out4,'y',ylen,yid)
status = NF_DEF_DIM(out4,'z',NSOLD,zid)
status = NF_DEF_DIM(out4,'tstep',NF_unlimited,tstepid)
outdim2(1) = xid
outdim2(2) = yid
outdim2(3) = zid
outdim2(4) = tstepid
outdim(1) = xid
outdim(2) = yid
outdim(3) = tstepid
status = NF_PUT_ATT_TEXT(out4,NF_GLOBAL,'Conventions',7,
& 'GDT 1.2')
status = NF_PUT_ATT_TEXT(out4,NF_GLOBAL,'file_name',33,
& 'noah.rerun.sub.torne.hrly.19'//YEAR_CH2//'.nc')
status = NF_PUT_ATT_TEXT(out4,NF_GLOBAL,'production',30,
& 'NOAH model run, DL on 20010312')
status = NF_DEF_VAR(out4,'lon',NF_FLOAT,1,xid,lonid4)
status = NF_DEF_VAR(out4,'lat',NF_FLOAT,1,yid,latid4)
status = NF_DEF_VAR(out4,'layerdepth',NF_FLOAT,1,zid,layerid)
status = NF_DEF_VAR(out4,'time',NF_FLOAT,1,tstepid,timeid4)
status = NF_DEF_VAR(out4,'timestp',NF_INT,1,tstepid,timestpid4)
status = NF_PUT_ATT_TEXT(out4,lonid4,'units',12,'degrees_east')
status = NF_PUT_ATT_REAL(out4,lonid4,'valid_min',NF_FLOAT,1,
& 18.125)
status = NF_PUT_ATT_REAL(out4,lonid4,'valid_max',NF_FLOAT,1,
& 25.125)
status = NF_PUT_ATT_TEXT(out4,lonid4,'long_name',9,'Longitude')
status = NF_PUT_ATT_TEXT(out4,latid4,'units',13,
& 'degrees_north')
status = NF_PUT_ATT_REAL(out4,latid4,'valid_min',NF_FLOAT,1,
& 65.875)
status = NF_PUT_ATT_REAL(out4,latid4,'valid_max',NF_FLOAT,1,
& 69.125)
status = NF_PUT_ATT_TEXT(out4,latid4,'long_name',8,'Latitude')
status = NF_PUT_ATT_TEXT(out4,layerid,'units',1,'m')
status = NF_PUT_ATT_TEXT(out4,layerid,'long_name',39,
& 'depth from surface to end of soil layer')
status = NF_PUT_ATT_TEXT(out4,timeid4,'units',33,
& 'seconds since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out4,timeid4,'calendar',9,
& 'gregorian')
status = NF_PUT_ATT_TEXT(out4,timeid4,'title',4,
& 'Time')
status = NF_PUT_ATT_TEXT(out4,timeid4,'long_name',9,
& 'Time axis')
status = NF_PUT_ATT_TEXT(out4,timeid4,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_PUT_ATT_TEXT(out4,timestpid4,'units',35,
& 'timesteps since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out4,timestpid4,'title',10,
& 'Time steps')
status = NF_PUT_ATT_REAL(out4,timestpid4,'tstep_sec',
& NF_FLOAT,1,3600.0)
status = NF_PUT_ATT_TEXT(out4,timestpid4,'long_name',14,
& 'Time step axis')
status = NF_PUT_ATT_TEXT(out4,timestpid4,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_DEF_VAR(out4,'SoilMoist',NF_FLOAT,4,outdim2,
& SoilMoistid)
status = NF_PUT_ATT_TEXT(out4,SoilMoistid,'units',6,'kg/m^2')
status = NF_PUT_ATT_TEXT(out4,SoilMoistid,'axis',12,
& 'Fortran XYZT')
status = NF_PUT_ATT_REAL(out4,SoilMoistid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out4,SoilMoistid,'associate',14,
& 'lon lat z time')
status = NF_DEF_VAR(out4,'LSoilMoist',NF_FLOAT,4,outdim2,
& LSoilMoistid)
status = NF_PUT_ATT_TEXT(out4,LSoilMoistid,'units',6,'kg/m^2')
status = NF_PUT_ATT_TEXT(out4,LSoilMoistid,'axis',12,
& 'Fortran XYZT')
status = NF_PUT_ATT_REAL(out4,LSoilMoistid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out4,LSoilMoistid,'associate',14,
& 'lon lat z time')
status = NF_DEF_VAR(out4,'SoilTemp',NF_FLOAT,4,outdim2,
& SoilTempid)
status = NF_PUT_ATT_TEXT(out4,SoilTempid,'units',1,'K')
status = NF_PUT_ATT_TEXT(out4,SoilTempid,'axis',12,
& 'Fortran XYZT')
status = NF_PUT_ATT_REAL(out4,SoilTempid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out4,SoilTempid,'associate',14,
& 'lon lat z time')
status = NF_DEF_VAR(out4,'SoilWet',NF_FLOAT,3,outdim,
& SoilWetid)
status = NF_PUT_ATT_TEXT(out4,SoilWetid,'units',1,'-')
status = NF_PUT_ATT_TEXT(out4,SoilWetid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out4,SoilWetid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out4,SoilWetid,'associate',12,
& 'lon lat time')
status = NF_ENDDEF(out4)
out5name = 'noah.rerun.eva.torne.hrly.19'//YEAR_CH2//'.nc'
flag = 0
status = NF_CREATE(out5name,flag,out5)
status = NF_DEF_DIM(out5,'x',xlen,xid)
status = NF_DEF_DIM(out5,'y',ylen,yid)
status = NF_DEF_DIM(out5,'tstep',NF_unlimited,tstepid)
outdim(1) = xid
outdim(2) = yid
outdim(3) = tstepid
status = NF_PUT_ATT_TEXT(out5,NF_GLOBAL,'Conventions',7,
& 'GDT 1.2')
status = NF_PUT_ATT_TEXT(out5,NF_GLOBAL,'file_name',33,
& 'noah.rerun.eva.torne.hrly.19'//YEAR_CH2//'.nc')
status = NF_PUT_ATT_TEXT(out5,NF_GLOBAL,'production',30,
& 'NOAH model run, DL on 20010312')
status = NF_DEF_VAR(out5,'lon',NF_FLOAT,1,xid,lonid5)
status = NF_DEF_VAR(out5,'lat',NF_FLOAT,1,yid,latid5)
status = NF_DEF_VAR(out5,'time',NF_FLOAT,1,tstepid,timeid5)
status = NF_DEF_VAR(out5,'timestp',NF_INT,1,tstepid,timestpid5)
status = NF_PUT_ATT_TEXT(out5,lonid5,'units',12,'degrees_east')
status = NF_PUT_ATT_REAL(out5,lonid5,'valid_min',NF_FLOAT,1,
& 18.125)
status = NF_PUT_ATT_REAL(out5,lonid5,'valid_max',NF_FLOAT,1,
& 25.125)
status = NF_PUT_ATT_TEXT(out5,lonid5,'long_name',9,'Longitude')
status = NF_PUT_ATT_TEXT(out5,latid5,'units',13,
& 'degrees_north')
status = NF_PUT_ATT_REAL(out5,latid5,'valid_min',NF_FLOAT,1,
& 65.875)
status = NF_PUT_ATT_REAL(out5,latid5,'valid_max',NF_FLOAT,1,
& 69.125)
status = NF_PUT_ATT_TEXT(out5,latid5,'long_name',8,'Latitude')
status = NF_PUT_ATT_TEXT(out5,timeid5,'units',33,
& 'seconds since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out5,timeid5,'calendar',9,
& 'gregorian')
status = NF_PUT_ATT_TEXT(out5,timeid5,'title',4,
& 'Time')
status = NF_PUT_ATT_TEXT(out5,timeid5,'long_name',9,
& 'Time axis')
status = NF_PUT_ATT_TEXT(out5,timeid5,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_PUT_ATT_TEXT(out5,timestpid5,'units',35,
& 'timesteps since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out5,timestpid5,'title',10,
& 'Time steps')
status = NF_PUT_ATT_REAL(out5,timestpid5,'tstep_sec',
& NF_FLOAT,1,3600.0)
status = NF_PUT_ATT_TEXT(out5,timestpid5,'long_name',14,
& 'Time step axis')
status = NF_PUT_ATT_TEXT(out5,timestpid5,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_DEF_VAR(out5,'ECanop',NF_FLOAT,3,outdim,
& ECanopid)
status = NF_PUT_ATT_TEXT(out5,ECanopid,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out5,ECanopid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out5,ECanopid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out5,ECanopid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out5,'ESoil',NF_FLOAT,3,outdim,
& ESoilid)
status = NF_PUT_ATT_TEXT(out5,ESoilid,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out5,ESoilid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out5,ESoilid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out5,ESoilid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out5,'TVeg',NF_FLOAT,3,outdim,
& TVegid)
status = NF_PUT_ATT_TEXT(out5,TVegid,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out5,TVegid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out5,TVegid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out5,TVegid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out5,'RootMoist',NF_FLOAT,3,outdim,
& RootMoistid)
status = NF_PUT_ATT_TEXT(out5,RootMoistid,'units',6,'kg/m^2')
status = NF_PUT_ATT_TEXT(out5,RootMoistid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out5,RootMoistid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out5,RootMoistid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out5,'CanopInt',NF_FLOAT,3,outdim,
& CanopIntid)
status = NF_PUT_ATT_TEXT(out5,CanopIntid,'units',6,'kg/m^2')
status = NF_PUT_ATT_TEXT(out5,CanopIntid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out5,CanopIntid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out5,CanopIntid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out5,'SubSnow',NF_FLOAT,3,outdim,
& SubSnowid)
status = NF_PUT_ATT_TEXT(out5,SubSnowid,'units',7,'kg/m^2s')
status = NF_PUT_ATT_TEXT(out5,SubSnowid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out5,SubSnowid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out5,SubSnowid,'associate',12,
& 'lon lat time')
status = NF_ENDDEF(out5)
out6name = 'noah.rerun.csp.torne.hrly.19'//YEAR_CH2//'.nc'
flag = 0
status = NF_CREATE(out6name,flag,out6)
status = NF_DEF_DIM(out6,'x',xlen,xid)
status = NF_DEF_DIM(out6,'y',ylen,yid)
status = NF_DEF_DIM(out6,'tstep',NF_unlimited,tstepid)
outdim(1) = xid
outdim(2) = yid
outdim(3) = tstepid
status = NF_PUT_ATT_TEXT(out6,NF_GLOBAL,'Conventions',7,
& 'GDT 1.2')
status = NF_PUT_ATT_TEXT(out6,NF_GLOBAL,'file_name',33,
& 'noah.rerun.csp.torne.hrly.19'//YEAR_CH2//'.nc')
status = NF_PUT_ATT_TEXT(out6,NF_GLOBAL,'production',30,
& 'NOAH model run, DL on 20010312')
status = NF_DEF_VAR(out6,'lon',NF_FLOAT,1,xid,lonid6)
status = NF_DEF_VAR(out6,'lat',NF_FLOAT,1,yid,latid6)
status = NF_DEF_VAR(out6,'time',NF_FLOAT,1,tstepid,timeid6)
status = NF_DEF_VAR(out6,'timestp',NF_INT,1,tstepid,timestpid6)
status = NF_PUT_ATT_TEXT(out6,lonid6,'units',12,'degrees_east')
status = NF_PUT_ATT_REAL(out6,lonid6,'valid_min',NF_FLOAT,1,
& 18.125)
status = NF_PUT_ATT_REAL(out6,lonid6,'valid_max',NF_FLOAT,1,
& 25.125)
status = NF_PUT_ATT_TEXT(out6,lonid6,'long_name',9,'Longitude')
status = NF_PUT_ATT_TEXT(out6,latid6,'units',13,
& 'degrees_north')
status = NF_PUT_ATT_REAL(out6,latid6,'valid_min',NF_FLOAT,1,
& 65.875)
status = NF_PUT_ATT_REAL(out6,latid6,'valid_max',NF_FLOAT,1,
& 69.125)
status = NF_PUT_ATT_TEXT(out6,latid6,'long_name',8,'Latitude')
status = NF_PUT_ATT_TEXT(out6,timeid6,'units',33,
& 'seconds since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out6,timeid6,'calendar',9,
& 'gregorian')
status = NF_PUT_ATT_TEXT(out6,timeid6,'title',4,
& 'Time')
status = NF_PUT_ATT_TEXT(out6,timeid6,'long_name',9,
& 'Time axis')
status = NF_PUT_ATT_TEXT(out6,timeid6,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_PUT_ATT_TEXT(out6,timestpid6,'units',35,
& 'timesteps since 1979-01-01 00:00:00')
status = NF_PUT_ATT_TEXT(out6,timestpid6,'title',10,
& 'Time steps')
status = NF_PUT_ATT_REAL(out6,timestpid6,'tstep_sec',
& NF_FLOAT,1,3600.0)
status = NF_PUT_ATT_TEXT(out6,timestpid6,'long_name',14,
& 'Time step axis')
status = NF_PUT_ATT_TEXT(out6,timestpid6,'time_origin',20,
& '1979-JAN-01 00:00:00')
status = NF_DEF_VAR(out6,'SnowFrac',NF_FLOAT,3,outdim,
& SnowFracid)
status = NF_PUT_ATT_TEXT(out6,SnowFracid,'units',1,'-')
status = NF_PUT_ATT_TEXT(out6,SnowFracid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out6,SnowFracid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out6,SnowFracid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out6,'SnowDepth',NF_FLOAT,3,outdim,
& SnowDepthid)
status = NF_PUT_ATT_TEXT(out6,SnowDepthid,'units',1,'m')
status = NF_PUT_ATT_TEXT(out6,SnowDepthid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out6,SnowDepthid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out6,SnowDepthid,'associate',12,
& 'lon lat time')
status = NF_DEF_VAR(out6,'SAlbedo',NF_FLOAT,3,outdim,
& SAlbedoid)
status = NF_PUT_ATT_TEXT(out6,SAlbedoid,'units',1,'-')
status = NF_PUT_ATT_TEXT(out6,SAlbedoid,'axis',11,
& 'Fortran XYT')
status = NF_PUT_ATT_REAL(out6,SAlbedoid,'missing_value',
& NF_FLOAT,1,1.0E+20)
status = NF_PUT_ATT_TEXT(out6,SAlbedoid,'associate',12,
& 'lon lat time')
status = NF_ENDDEF(out6)
DO NT = 1,tsteplen
timestp = timestp + 1
status = NF_PUT_VARA_REAL(out1,timeid1,NT,1,timestp*3600.0)
status = NF_PUT_VARA_INT(out1,timestpid1,NT,1,timestp)
status = NF_PUT_VARA_REAL(out2,timeid2,NT,1,timestp*3600.0)
status = NF_PUT_VARA_INT(out2,timestpid2,NT,1,timestp)
status = NF_PUT_VARA_REAL(out3,timeid3,NT,1,timestp*3600.0)
status = NF_PUT_VARA_INT(out3,timestpid3,NT,1,timestp)
status = NF_PUT_VARA_REAL(out4,timeid4,NT,1,timestp*3600.0)
status = NF_PUT_VARA_INT(out4,timestpid4,NT,1,timestp)
status = NF_PUT_VARA_REAL(out5,timeid5,NT,1,timestp*3600.0)
status = NF_PUT_VARA_INT(out5,timestpid5,NT,1,timestp)
status = NF_PUT_VARA_REAL(out6,timeid6,NT,1,timestp*3600.0)
status = NF_PUT_VARA_INT(out6,timestpid6,NT,1,timestp)
C CALCULATE JULIAN DAY, CALCULATE WEIGTHS AND INTERPOLATE FROM
C MONTHLY TO DAILY
JULDAY = mod(INT((NT-1)/24),365) + 1
CALL CALC_WEIGTHS(JULDAY,W1,W2,M1,M2)
DO J = 1,NY
DO I = 1,NX
ALBEDO_D(I,J) = W1 * ALBEDO(M1,I,J) +
& W2 * ALBEDO(M2,I,J)
SHDFAC_D(I,J) = W1 * SHDFAC(M1,I,J) +
& W2 * SHDFAC(M2,I,J)
END DO
END DO
start(1) = 1
start(2) = 1
IF (NT .LT. tsteplen) THEN
start(3) = NT+1
ELSE
start(3) = NT
END IF
count(1) = xlen
count(2) = ylen
count(3) = 1
status = NF_INQ_VARID(ncid,'SWdown',SWdownid)
status = NF_INQ_VAR(ncid,SWdownid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,SWdownid,start,count,
& SWdownplus)
status = NF_INQ_VARID(ncid,'Rainf',Rainfid)
status = NF_INQ_VAR(ncid,Rainfid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Rainfid,start,count,
& Rainfplus)
status = NF_INQ_VARID(ncid,'Snowf',Snowfid)
status = NF_INQ_VAR(ncid,Snowfid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Snowfid,start,count,
& Snowfplus)
status = NF_INQ_VARID(ncid,'LWdown',LWdownid)
status = NF_INQ_VAR(ncid,LWdownid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,LWdownid,start,count,
& LWdownplus)
status = NF_INQ_VARID(ncid,'PSurf',PSurfid)
status = NF_INQ_VAR(ncid,PSurfid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,PSurfid,start,count,
& PSurfplus)
status = NF_INQ_VARID(ncid,'Tair',Tairid)
status = NF_INQ_VAR(ncid,Tairid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Tairid,start,count,
& Tairplus)
status = NF_INQ_VARID(ncid,'Qair',Qairid)
status = NF_INQ_VAR(ncid,Qairid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Qairid,start,count,
& Qairplus)
status = NF_INQ_VARID(ncid,'Wind_N',Wind_Nid)
status = NF_INQ_VAR(ncid,Wind_Nid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Wind_Nid,start,count,
& Wind_Nplus)
status = NF_INQ_VARID(ncid,'Wind_E',Wind_Eid)
status = NF_INQ_VAR(ncid,Wind_Eid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Wind_Eid,start,count,
& Wind_Eplus)
IF (FIRST_STEP) THEN
start(1) = 1
start(2) = 1
start(3) = 1
count(1) = xlen
count(2) = ylen
count(3) = 1
status = NF_INQ_VARID(ncid,'SWdown',SWdownid)
status = NF_INQ_VAR(ncid,SWdownid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,SWdownid,start,count,
& SWdownnow)
status = NF_INQ_VARID(ncid,'Rainf',Rainfid)
status = NF_INQ_VAR(ncid,Rainfid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Rainfid,start,count,
& Rainfnow)
status = NF_INQ_VARID(ncid,'Snowf',Snowfid)
status = NF_INQ_VAR(ncid,Snowfid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Snowfid,start,count,
& Snowfnow)
status = NF_INQ_VARID(ncid,'LWdown',LWdownid)
status = NF_INQ_VAR(ncid,LWdownid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,LWdownid,start,count,
& LWdownnow)
status = NF_INQ_VARID(ncid,'PSurf',PSurfid)
status = NF_INQ_VAR(ncid,PSurfid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,PSurfid,start,count,
& PSurfnow)
status = NF_INQ_VARID(ncid,'Tair',Tairid)
status = NF_INQ_VAR(ncid,Tairid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Tairid,start,count,
& Tairnow)
status = NF_INQ_VARID(ncid,'Qair',Qairid)
status = NF_INQ_VAR(ncid,Qairid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Qairid,start,count,
& Qairnow)
status = NF_INQ_VARID(ncid,'Wind_N',Wind_Nid)
status = NF_INQ_VAR(ncid,Wind_Nid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Wind_Nid,start,count,
& Wind_Nnow)
status = NF_INQ_VARID(ncid,'Wind_E',Wind_Eid)
status = NF_INQ_VAR(ncid,Wind_Eid,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VARA_REAL(ncid,Wind_Eid,start,count,
& Wind_Enow)
status = NF_INQ_VARID(ncid,'lon',lon2id)
status = NF_INQ_VAR(ncid,lon2id,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VAR_REAL(ncid,lon2id,lon)
status = NF_PUT_VAR_REAL(out1,lonid1,lon)
status = NF_PUT_VAR_REAL(out2,lonid2,lon)
status = NF_PUT_VAR_REAL(out3,lonid3,lon)
status = NF_PUT_VAR_REAL(out4,lonid4,lon)
status = NF_PUT_VAR_REAL(out5,lonid5,lon)
status = NF_PUT_VAR_REAL(out6,lonid6,lon)
status = NF_INQ_VARID(ncid,'lat',lat2id)
status = NF_INQ_VAR(ncid,lat2id,name,xtype,ndims,
& dimids,natts)
status = NF_GET_VAR_REAL(ncid,lat2id,lat)
status = NF_PUT_VAR_REAL(out1,latid1,lat)
status = NF_PUT_VAR_REAL(out2,latid2,lat)
status = NF_PUT_VAR_REAL(out3,latid3,lat)
status = NF_PUT_VAR_REAL(out4,latid4,lat)
status = NF_PUT_VAR_REAL(out5,latid5,lat)
status = NF_PUT_VAR_REAL(out6,latid6,lat)
status = NF_PUT_VAR_REAL(out4,layerid,layerdepth)
SWdownminus = SWdownnow
LWdownminus = LWdownnow
PSurfminus = PSurfnow
Tairminus = Tairnow
Qairminus = Qairnow
WIND_Nminus = WIND_Nnow
Wind_Eminus = Wind_Enow
DO J = 1,NY
DO I = 1,NX
SOILM(I,J) = 0.0
IF (LAND_SEA(I,J) .EQ. 1) THEN
DO K = 1, NSOIL(I,J)
SOILM(I,J)=SOILM(I,J)+SMC(K,I,J)*
& SOILDEPTH(K,I,J)
END DO
END IF
END DO
END DO
FIRST_STEP = .FALSE.
END IF
HOUR_SOLDN = 0.0
HOUR_LWDN = 0.0
HOUR_LWDN_MODEL = 0.0
HOUR_ETA = 0.0
HOUR_H = 0.0
HOUR_S = 0.0
HOUR_RUNOFF1 = 0.0
HOUR_RUNOFF2 = 0.0
HOUR_SNMAX = 0.0
SoilMoist = 0.0
LSoilMoist = 0.0
SoilTemp = 0.0
SoilWet = 0.0
HOUR_ALBEDO = 0.0
HOUR_T = 0.0
SNEQV_ALMA = 0.0
SnowT = 0.0
HOUR_CMC = 0.0
HOUR_EDIR1 = 0.0
HOUR_EC1 = 0.0
HOUR_ETT1 = 0.0
HOUR_ESNOW = 0.0
RootMoist = 0.0
SnowFrac = 0.0
HOUR_SNOWDEPTH = 0.0
CMC_old = CMC
SOILM_old = SOILM
SNEQV_old = SNEQV
DO SUB_DT = 1,4
DO J = 1,NY
DO I = 1,NX
C DO J = 8,8
C DO I = 15,15
IF (LAND_SEA(I,J) .EQ. 1) THEN
DT_TAIR = (Tairminus(I,J)* W_minus(SUB_DT) +
& Tairnow(I,J) * W_now(SUB_DT) +
& Tairplus(I,J) * W_plus(SUB_DT)) *
& 4.0*Tairnow(I,J) / (0.5*Tairminus(I,J) +
& 3.0*Tairnow(I,J) + 0.5*Tairplus(I,J))
DT_SPFH = (Qairminus(I,J)* W_minus(SUB_DT) +
& Qairnow(I,J) * W_now(SUB_DT) +
& Qairplus(I,J) * W_plus(SUB_DT)) *
& 4.0*Qairnow(I,J) / (0.5*Qairminus(I,J) +
& 3.0*Qairnow(I,J) + 0.5*Qairplus(I,J))
DT_PSFC = (PSurfminus(I,J)* W_minus(SUB_DT) +
& PSurfnow(I,J) * W_now(SUB_DT) +
& PSurfplus(I,J) * W_plus(SUB_DT)) *
& 4.*PSurfnow(I,J) / (.5*PSurfminus(I,J) +
& 3.0*PSurfnow(I,J) + 0.5*PSurfplus(I,J))
IF ((0.5*Wind_Eminus(I,J)+3.0*Wind_Enow(I,J)+
& 0.5*Wind_Eplus(I,J)) .GT. 0.01) THEN
DT_UWIND=(Wind_Eminus(I,J)*W_minus(SUB_DT)+
& Wind_Enow(I,J) * W_now(SUB_DT) +
& Wind_Eplus(I,J) * W_plus(SUB_DT)) *
& 4.*Wind_Enow(I,J)/(.5*Wind_Eminus(I,J)+
& 3.0*Wind_Enow(I,J)+0.5*Wind_Eplus(I,J))
ELSE
DT_UWIND=0.0
END IF
IF ((0.5*Wind_Nminus(I,J)+3.0*Wind_Nnow(I,J) +
& 0.5*Wind_Nplus(I,J)) .GT. 0.01) THEN
DT_VWIND=(Wind_Nminus(I,J)*W_minus(SUB_DT)+
& Wind_Nnow(I,J) * W_now(SUB_DT) +
& Wind_Nplus(I,J) * W_plus(SUB_DT)) *
& 4.*Wind_Nnow(I,J)/(.5*Wind_Nminus(I,J)+
& 3.0*Wind_Nnow(I,J)+0.5*Wind_Nplus(I,J))
ELSE
DT_VWIND=0.0
END IF
DT_LWDN = (LWdownminus(I,J)* W_minus(SUB_DT) +
& LWdownnow(I,J) * W_now(SUB_DT) +
& LWdownplus(I,J) * W_plus(SUB_DT)) *
& 4.0*LWdownnow(I,J)/(0.5*LWdownminus(I,J)+
& 3.0*LWdownnow(I,J) + 0.5*LWdownplus(I,J))
IF ( (0.5*SWdownminus(I,J)+3.0*SWdownnow(I,J)+
& 0.5*SWdownplus(I,J)) .GT. 0.01) THEN
DT_SOLDN=(SWdownminus(I,J)*W_minus(SUB_DT)+
& SWdownnow(I,J) * W_now(SUB_DT) +
& SWdownplus(I,J) * W_plus(SUB_DT)) *
& 4.*SWdownnow(I,J)/(.5*SWdownminus(I,J)+
& 3.0*SWdownnow(I,J)+0.5*SWdownplus(I,J))
ELSE
DT_SOLDN = 0.0
END IF
DT_PRCP = (Rainfnow(I,J)+Snowfnow(I,J)) *
& W_PRECIP(SUB_DT)
C DRIVER STEP 6
C CALCULATE CH (EXCHANGE COEFFICIENT)
CALL QDATAP(DT_TAIR,DT_PSFC,DT_SPFH_SAT)
IF (DT_SPFH .GE. DT_SPFH_SAT) THEN
DT_SPFH = DT_SPFH_SAT
END IF
C CALCULATE SLOPE OF SAT SPECIFIC HUMIDITY CURVE FOR PENMAN: DQSDT2
DQSDT2 = DQSDT(DT_TAIR, DT_PSFC)
TH2 = DT_TAIR + 0.0098 * Z
SFCSPD = SQRT(DT_UWIND*DT_UWIND+
& DT_VWIND*DT_VWIND)
C DRIVER STEP 7
C CALL LAND-SURFACE PHYSICS
CALL SFLX(ICE,DT,Z,NSOIL(I,J),
I SOILDEPTH(1,I,J),DT_LWDN,DT_SOLDN,
I DT_PSFC,DT_PRCP,DT_TAIR,TH2,DT_SPFH,
I DT_SPFH_SAT,DQSDT2,VEGTYP(I,J),SOILTYP(I,J),
I SLOPETYP(I,J),SHDFAC_D(I,J),PTU,
I TBOT(I,J),ALBEDO_D(I,J),SNOALB,SFCSPD,
2 CMC(I,J),T1(I,J),STC(1,I,J),SMC(1,I,J),
I SH2O(1,I,J),SNOWH(I,J),SNEQV(I,J),CH(I,J),
O CM(I,J),ETP(I,J),ETA(I,J),H(I,J),S(I,J),
O RUNOFF1(I,J),RUNOFF2(I,J),Q1(I,J),SNMAX(I,J),
O ALBSNO(I,J),SOILW(I,J),SOILM(I,J),
O SMCWLT,SMCDRY,SMCREF,SMCMAX,EDIR1(I,J),
O EC1(I,J),ETT1(I,J),ESNOW(I,J),RM,SCOVER)
C Energy balance output variables
HOUR_SOLDN(I,J)= HOUR_SOLDN(I,J)+(1.0-
& ALBSNO(I,J))* DT_SOLDN/4.0
HOUR_LWDN(I,J) = HOUR_LWDN(I,J)+DT_LWDN/4.0
HOUR_LWDN_MODEL(I,J) = HOUR_LWDN_MODEL(I,J) +
& 5.67E-8 * T1(I,J)**4.0 / 4.0
HOUR_ETA(I,J) = HOUR_ETA(I,J) - ETA(I,J)/4.0
HOUR_H(I,J) = HOUR_H(I,J) - H(I,J)/4.0
HOUR_S(I,J) = HOUR_S(I,J) + S(I,J)/4.0
C Water balance output variables
HOUR_RUNOFF1(I,J) = HOUR_RUNOFF1(I,J) -
& RUNOFF1(I,J) * 1000.0 / 4.0
HOUR_RUNOFF2(I,J) = HOUR_RUNOFF2(I,J) -
& RUNOFF2(I,J) * 1000.0 / 4.0
HOUR_SNMAX(I,J) = HOUR_SNMAX(I,J) +
& SNMAX(I,J)
C Surface State Variables
DO K = 1, NSOLD
SoilMoist(I,J,K) = SoilMoist(K,I,J) +
& 1000.0*SMC(K,I,J)*
& SOILDEPTH(K,I,J)/4.0
LSoilMoist(I,J,K) = LSoilMoist(K,I,J)+
& 1000.0*SH2O(K,I,J)*
& SOILDEPTH(K,I,J)/4.0
SoilTemp(I,J,K) = SoilTemp(K,I,J) +
& STC(K,I,J) / 4.0
END DO
SoilWet(I,J) = SoilWet(I,J) +
& SOILW(I,J) / 4.0
HOUR_ALBEDO(I,J) = HOUR_ALBEDO(I,J) +
& ALBSNO(I,J) / 4.0
HOUR_T(I,J) = HOUR_T(I,J) + T1(I,J) / 4.0
SNEQV_ALMA(I,J) = SNEQV_ALMA(I,J) +
& SNEQV(I,J)*1000.0/4.0
SnowT(I,J) = SNOWT(I,J) +
& 0.5*(T1(I,J)+STC(1,I,J))/4.0
C Evaporation Components
HOUR_CMC(I,J) = CMC(I,J) / 4.0
HOUR_EDIR1(I,J) = HOUR_EDIR1(I,J) -
& 1000.0 * EDIR1(I,J)/4.0
HOUR_EC1(I,J) = HOUR_EC1(I,J) -
& 1000.0 * EC1(I,J)/4.0
HOUR_ETT1(I,J) = HOUR_ETT1(I,J) -
& 1000.0 * ETT1(I,J)/4.0
HOUR_ESNOW(I,J) = HOUR_ESNOW(I,J) -
& ESNOW(I,J) / 4.0
RootMoist(I,J) = RootMoist(I,J) + RM/4.0
C Cold Season Processes
SnowFrac(I,J) = SnowFrac(I,J)+SCOVER/4.0
HOUR_SNOWDEPTH(I,J)=HOUR_SNOWDEPTH(I,J)+
& SNOWH(I,J) / 4.0
END IF
END DO
END DO
END DO
start(1) = 1
start(2) = 1
start(3) = NT
count(1) = xlen
count(2) = ylen
count(3) = 1
start2(1) = 1
start2(2) = 1
start2(3) = 1
start2(4) = NT
count2(1) = xlen
count2(2) = ylen
count2(3) = NSOLD
count2(4) = 1
HOUR_LWDN_BUDGET = HOUR_LWDN - HOUR_LWDN_MODEL
WHERE (LAND_SEA .EQ. 0)
HOUR_SOLDN = 1.0E+20
HOUR_LWDN_BUDGET = 1.0E+20
HOUR_ETA = 1.0E+20
HOUR_H = 1.0E+20
HOUR_S = 1.0E+20
END WHERE
status = NF_PUT_VARA_REAL(out1,SWnetid,start,count,
& HOUR_SOLDN)
status = NF_PUT_VARA_REAL(out1,LWnetid,start,count,
& HOUR_LWDN_BUDGET)
status = NF_PUT_VARA_REAL(out1,Qleid,start,count,
& HOUR_ETA)
status = NF_PUT_VARA_REAL(out1,Qhid,start,count,HOUR_H)
status = NF_PUT_VARA_REAL(out1,Qgid,start,count,HOUR_S)
HOUR_EVAP = HOUR_ETA / 2.501000E+6
DelSoil = 1000.0 * (SOILM - SOILM_old)
DelInter = 1000.0 * (CMC - CMC_old)
DelSWE = 1000.0 * (SNEQV - SNEQV_old)
HOUR_SNMAX = 1000.0 * HOUR_SNMAX / (4.0 * DT)
WHERE (LAND_SEA .EQ. 0)
Snowfnow = 1.0E+20
Rainfnow = 1.0E+20
HOUR_EVAP = 1.0E+20
HOUR_RUNOFF1 = 1.0E+20
HOUR_RUNOFF2 = 1.0E+20
HOUR_SNMAX = 1.0E+20
DelSoil = 1.0E+20
DelSWE = 1.0E+20
DelInter = 1.0E+20
END WHERE
status = NF_PUT_VARA_REAL(out2,Snowf2id,start,count,
& Snowfnow)
status = NF_PUT_VARA_REAL(out2,Rainf2id,start,count,
& Rainfnow)
status = NF_PUT_VARA_REAL(out2,Evapid,start,count,
& HOUR_EVAP)
status = NF_PUT_VARA_REAL(out2,Qsid,start,count,
& HOUR_RUNOFF1)
status = NF_PUT_VARA_REAL(out2,Qsbid,start,count,
& HOUR_RUNOFF2)
status = NF_PUT_VARA_REAL(out2,Qsmid,start,count,
& HOUR_SNMAX)
status = NF_PUT_VARA_REAL(out2,DelSoilid,start,count,
& DelSoil)
status = NF_PUT_VARA_REAL(out2,DelSWEid,start,count,
& DelSWE)
status = NF_PUT_VARA_REAL(out2,DelInterid,start,count,
& DelInter)
WHERE (SNEQV_ALMA .GT. 0.0)
SnowT = min(273.15,SnowT)
END WHERE
WHERE (SNEQV_ALMA .EQ. 0.0)
SnowT = 1.0E+20
END WHERE
WHERE (LAND_SEA .EQ. 0)
SnowT = 1.0E+20
HOUR_T = 1.0E+20
SNEQV_ALMA = 1.0E+20
HOUR_ALBEDO = 1.0E+20
END WHERE
status = NF_PUT_VARA_REAL(out3,SnowTid,start,count,
& SnowT)
status = NF_PUT_VARA_REAL(out3,BaresoilTid,start,count,
& HOUR_T)
status = NF_PUT_VARA_REAL(out3,RadTid,start,count,
& HOUR_T)
status = NF_PUT_VARA_REAL(out3,VegTid,start,count,
& HOUR_T)
status = NF_PUT_VARA_REAL(out3,AvgSurfTid,start,count,
& HOUR_T)
status = NF_PUT_VARA_REAL(out3,SWEid,start,count,
& SNEQV_ALMA)
status = NF_PUT_VARA_REAL(out3,Albedoid,start,count,
& HOUR_ALBEDO)
DO J = 1, NY
DO I = 1, NX
IF (LAND_SEA(I,J) .EQ. 0) THEN
DO K = 1, NSOLD
SoilMoist(I,J,K) = 1.0E+20
LSoilMoist(I,J,K) = 1.0E+20
SoilTemp(I,J,K) = 1.0E+20
END DO
SoilWet(I,J) = 1.0E+20
END IF
END DO
END DO
status = NF_PUT_VARA_REAL(out4,SoilMoistid,start2,count2,
& SoilMoist)
status = NF_PUT_VARA_REAL(out4,LSoilMoistid,start2,count2,
& LSoilMoist)
status = NF_PUT_VARA_REAL(out4,SoilTempid,start2,count2,
& SoilTemp)
status = NF_PUT_VARA_REAL(out4,SoilWetid,start,count,
& SoilWet)
WHERE (LAND_SEA .EQ. 0)
HOUR_EC1 = 1.0E+20
HOUR_EDIR1 = 1.0E+20
HOUR_ETT1 = 1.0E+20
HOUR_CMC = 1.0E+20
HOUR_ESNOW = 1.0E+20
END WHERE
status = NF_PUT_VARA_REAL(out5,ECanopid,start,count,
& HOUR_EC1)
status = NF_PUT_VARA_REAL(out5,ESoilid,start,count,
& HOUR_EDIR1)
status = NF_PUT_VARA_REAL(out5,TVegid,start,count,
& HOUR_ETT1)
status = NF_PUT_VARA_REAL(out5,CanopIntid,start,count,
& HOUR_CMC)
status = NF_PUT_VARA_REAL(out5,SubSnowid,start,count,
& HOUR_ESNOW)
status = NF_PUT_VARA_REAL(out5,RootMoistid,start,count,
& RootMoist)
SAlbedo = 0.0
WHERE (SNEQV_ALMA .GT. 0.0)
SAlbedo = HOUR_ALBEDO + (1.0-SnowFrac)*
& (SNOALB-HOUR_ALBEDO)*(1.0-SHDFAC_D)
END WHERE
WHERE (LAND_SEA .EQ. 0)
SnowFrac = 1.0E+20
HOUR_SNOWDEPTH = 1.0E+20
SAlbedo = 1.0E+20
END WHERE
status = NF_PUT_VARA_REAL(out6,SnowFracid,start,count,
& SnowFrac)
status = NF_PUT_VARA_REAL(out6,SnowDepthid,start,count,
& HOUR_SNOWDEPTH)
status = NF_PUT_VARA_REAL(out6,SAlbedoid,start,count,
& SAlbedo)
SWdownminus = SWdownnow
LWdownminus = LWdownnow
PSurfminus = PSurfnow
Tairminus = Tairnow
Qairminus = Qairnow
WIND_Nminus = WIND_Nnow
Wind_Eminus = Wind_Enow
SWdownnow = SWdownplus
Rainfnow = Rainfplus
Snowfnow = Snowfplus
LWdownnow = LWdownplus
PSurfnow = PSurfplus
Tairnow = Tairplus
Qairnow = Qairplus
WIND_Nnow = WIND_Nplus
Wind_Enow = Wind_Eplus
END DO
status = NF_CLOSE(out1)
status = NF_CLOSE(out2)
status = NF_CLOSE(out3)
status = NF_CLOSE(out4)
status = NF_CLOSE(out5)
status = NF_CLOSE(out6)
status = NF_CLOSE(ncid)
WRITE(*,*) 'writing restart files'
OPEN(99, FILE = 'restart.T1'//YEAR_CH2, FORM = 'UNFORMATTED')
WRITE(99) ((T1(I,J),I=1,NX),J=1,NY)
CLOSE(99)
OPEN(99, FILE = 'restart.STC'//YEAR_CH2, FORM = 'UNFORMATTED')
WRITE(99) (((STC(K,I,J),K=1,NSOLD),I=1,NX),J=1,NY)
CLOSE(99)
OPEN(99, FILE = 'restart.SMC'//YEAR_CH2, FORM = 'UNFORMATTED')
WRITE(99) (((SMC(K,I,J),K=1,NSOLD),I=1,NX),J=1,NY)
CLOSE(99)
OPEN(99, FILE = 'restart.SH2O'//YEAR_CH2, FORM = 'UNFORMATTED')
WRITE(99) (((SH2O(K,I,J),K=1,NSOLD),I=1,NX),J=1,NY)
CLOSE(99)
OPEN(99, FILE = 'restart.CMC'//YEAR_CH2, FORM = 'UNFORMATTED')
WRITE(99) ((CMC(I,J),I=1,NX),J=1,NY)
CLOSE(99)
OPEN(99, FILE = 'restart.CH'//YEAR_CH2, FORM = 'UNFORMATTED')
WRITE(99) ((CH(I,J),I=1,NX),J=1,NY)
CLOSE(99)
OPEN(99, FILE = 'restart.CM'//YEAR_CH2, FORM = 'UNFORMATTED')
WRITE(99) ((CM(I,J),I=1,NX),J=1,NY)
CLOSE(99)
OPEN(99, FILE='restart.SNOWH'//YEAR_CH2,FORM ='UNFORMATTED')
WRITE(99) ((SNOWH(I,J),I=1,NX),J=1,NY)
CLOSE(99)
OPEN(99, FILE='restart.SNEQV'//YEAR_CH2,FORM='UNFORMATTED')
WRITE(99) ((SNEQV(I,J),I=1,NX),J=1,NY)
CLOSE(99)
END DO
END