MODULE AtmosModel USE Constants, ONLY: & r8, & r4, & i8, & pi USE Utils, ONLY :julday,lati USE Options, ONLY :lon_site, lat_site,icn_data,DELTAOUT !PK USE ReadNETCDF , ONLY : RunReadNETCDF,ncid_prec,ncid_temp,ncid_rhum,ncid_wind,ncid_swrd,ncid_lwrd,ncid_umes,ncid_pres IMPLICIT NONE PRIVATE REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gu (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gv (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gt (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gq (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: gps (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: dcupr_sib (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: dlspr_sib (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: radvbc_sib (: ,:)!calc em funcao de swdown REAL(KIND=r8), PUBLIC, ALLOCATABLE :: radnbc_sib (: ,:)!calc em funcao de swdown REAL(KIND=r8), PUBLIC, ALLOCATABLE :: radvdc_sib (: ,:)!calc em funcao de swdown REAL(KIND=r8), PUBLIC, ALLOCATABLE :: radndc_sib (: ,:)!calc em funcao de swdown REAL(KIND=r8), PUBLIC, ALLOCATABLE :: swdown(:,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: dlwbot_sib (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: cosz_sib (: ,:)!calc em funcao de latitude REAL(KIND=r8), PUBLIC, ALLOCATABLE :: ws (:,:,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: psb_g (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: bps_g (: ,:) REAL(KIND=r8), PUBLIC :: sigki REAL(KIND=r8), PUBLIC :: delsig REAL(KIND=r8), PUBLIC, ALLOCATABLE :: mskant8 (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: tseam (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: tsea (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: slrad (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: tsurf (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: qsurf (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: zorl (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: tmtx (:,:,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: qmtx (:,:,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: umtx (:,:,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: slhf (: ,:) REAL(KIND=r8), PUBLIC, ALLOCATABLE :: sshf (: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: latic(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: longc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: t02mc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: q02mc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: uvesc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: vvesc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: pslcc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: precc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: ocisc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: slhfc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: sshfc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: olisc(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: latip(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: longp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: t02mp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: q02mp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: uvesp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: vvesp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: pslcp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: precp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: ocisp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: slhfp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: sshfp(: ,:) REAL(KIND=r4), PUBLIC, ALLOCATABLE :: olisp(: ,:) INTEGER :: monday(12) INTEGER , PUBLIC :: monl(12)=(/31,28,31,30,31,30,31,31,30,31,30,31/) REAL (KIND=r8), PARAMETER :: pie = 3.1415926e0_r8! constant pi=3.1415926e0 REAL (KIND=r8), PARAMETER :: pai12 = pie/12.0_r8 ! pi/12 REAL (KIND=r8), PARAMETER :: pai = 3.14159265358979_r8! constant pi=3.1415926e0 LOGICAL :: OPENFILE=.TRUE. INTEGER :: ibMax INTEGER :: jbMax INTEGER :: nband INTEGER :: step PUBLIC :: InitAtmosModel PUBLIC :: RunAtmosModel2D PUBLIC :: FinalizeAtmosModel PUBLIC :: RunAtmosModel1D CONTAINS SUBROUTINE InitAtmosModel(ibMax_in,jbMax_in,nband_in) IMPLICIT NONE INTEGER, INTENT(IN ) :: ibMax_in INTEGER, INTENT(IN ) :: jbMax_in INTEGER, INTENT(IN ) :: nband_in ibMax = ibMax_in jbMax = jbMax_in nband = nband_in ALLOCATE(gu (1:ibMax ,jbMax) ) ALLOCATE(gv (1:ibMax ,jbMax) ) ALLOCATE(gt (1:ibMax ,jbMax) ) ALLOCATE(gq (1:ibMax ,jbMax) ) ALLOCATE(gps (1:ibMax ,jbMax) ) ALLOCATE(dcupr_sib (1:ibMax ,jbMax) ) ALLOCATE(dlspr_sib (1:ibMax ,jbMax) ) ALLOCATE(radvbc_sib (1:ibMax ,jbMax) ) ALLOCATE(radnbc_sib (1:ibMax ,jbMax) ) ALLOCATE(radvdc_sib (1:ibMax ,jbMax) ) ALLOCATE(radndc_sib (1:ibMax ,jbMax) ) ALLOCATE(dlwbot_sib (1:ibMax ,jbMax) ) ALLOCATE(swdown (1:ibMax ,jbMax) ) ALLOCATE(cosz_sib (1:ibMax ,jbMax) ) ALLOCATE(ws (1:ibMax,1:3,jbMax) ) ALLOCATE(psb_g (1:ibMax ,jbMax) ) ALLOCATE(bps_g (1:ibMax ,jbMax) ) ALLOCATE(mskant8 (1:ibMax ,jbMax) ) ALLOCATE(tseam (1:ibMax ,jbMax) ) ALLOCATE(tsea (1:ibMax ,jbMax) ) ALLOCATE(slrad (1:ibMax ,jbMax) ) ALLOCATE( slhf (1:ibMax ,jbMax) ) ALLOCATE( sshf (1:ibMax ,jbMax) ) ALLOCATE(tsurf (1:ibMax ,jbMax) ) ALLOCATE(qsurf (1:ibMax ,jbMax) ) ALLOCATE(zorl (1:ibMax ,jbMax) ) ALLOCATE(tmtx (1:ibMax,1:3,jbMax) ) ALLOCATE(qmtx (1:ibMax,1:3,jbMax) ) ALLOCATE(umtx (1:ibMax,1:4,jbMax) ) ALLOCATE(latic(1:ibMax ,jbMax) ) ALLOCATE(longc(1:ibMax ,jbMax) ) ALLOCATE(t02mc(1:ibMax ,jbMax) ) ALLOCATE(q02mc(1:ibMax ,jbMax) ) ALLOCATE(uvesc(1:ibMax ,jbMax) ) ALLOCATE(vvesc(1:ibMax ,jbMax) ) ALLOCATE(pslcc(1:ibMax ,jbMax) ) ALLOCATE(precc(1:ibMax ,jbMax) ) ALLOCATE(ocisc(1:ibMax ,jbMax) ) ALLOCATE(slhfc(1:ibMax ,jbMax) ) ALLOCATE(sshfc(1:ibMax ,jbMax) ) ALLOCATE(olisc(1:ibMax ,jbMax) ) ALLOCATE(latip(1:ibMax ,jbMax) ) ALLOCATE(longp(1:ibMax ,jbMax) ) ALLOCATE(t02mp(1:ibMax ,jbMax) ) ALLOCATE(q02mp(1:ibMax ,jbMax) ) ALLOCATE(uvesp(1:ibMax ,jbMax) ) ALLOCATE(vvesp(1:ibMax ,jbMax) ) ALLOCATE(pslcp(1:ibMax ,jbMax) ) ALLOCATE(precp(1:ibMax ,jbMax) ) ALLOCATE(ocisp(1:ibMax ,jbMax) ) ALLOCATE(slhfp(1:ibMax ,jbMax) ) ALLOCATE(sshfp(1:ibMax ,jbMax) ) ALLOCATE(olisp(1:ibMax ,jbMax) ) CALL InitRadtim(monl) END SUBROUTINE InitAtmosModel SUBROUTINE FinalizeAtmosModel() IMPLICIT NONE DEALLOCATE(gu ) DEALLOCATE(gv ) DEALLOCATE(gt ) DEALLOCATE(gq ) DEALLOCATE(gps ) DEALLOCATE(dcupr_sib ) DEALLOCATE(dlspr_sib ) DEALLOCATE(radvbc_sib) DEALLOCATE(radnbc_sib) DEALLOCATE(radvdc_sib) DEALLOCATE(radndc_sib) DEALLOCATE(dlwbot_sib) DEALLOCATE( swdown) DEALLOCATE(cosz_sib ) DEALLOCATE(ws ) DEALLOCATE(psb_g ) DEALLOCATE(bps_g ) DEALLOCATE(mskant8 ) DEALLOCATE(tseam ) DEALLOCATE(tsea ) DEALLOCATE(slrad ) DEALLOCATE( slhf ) DEALLOCATE( sshf ) DEALLOCATE(tsurf ) DEALLOCATE(qsurf ) DEALLOCATE(zorl ) DEALLOCATE(tmtx ) DEALLOCATE(qmtx ) DEALLOCATE(umtx ) DEALLOCATE(latic ) DEALLOCATE(longc ) DEALLOCATE(t02mc ) DEALLOCATE(q02mc ) DEALLOCATE(uvesc ) DEALLOCATE(vvesc ) DEALLOCATE(pslcc ) DEALLOCATE(precc ) DEALLOCATE(ocisc ) DEALLOCATE(slhfc ) DEALLOCATE(sshfc ) DEALLOCATE(olisc ) DEALLOCATE(latip ) DEALLOCATE(longp ) DEALLOCATE(t02mp ) DEALLOCATE(q02mp ) DEALLOCATE(uvesp ) DEALLOCATE(vvesp ) DEALLOCATE(pslcp ) DEALLOCATE(precp ) DEALLOCATE(ocisp ) DEALLOCATE(slhfp ) DEALLOCATE(sshfp ) DEALLOCATE(olisp ) END SUBROUTINE FinalizeAtmosModel SUBROUTINE RunAtmosModel1D(jdt,tod,nband,dt,idate,idatec,idatep,DELTAIN,start) IMPLICIT NONE REAL(KIND=r8), INTENT(IN ) :: tod INTEGER , INTENT(IN ) :: jdt INTEGER , INTENT(IN ) :: nband ! INTENT(IN ) :: nband REAL(KIND=r8), INTENT(IN ) :: dt ! INTENT(IN ) :: dt INTEGER , INTENT(IN ) :: idate (4) ! INTENT(IN ) :: idatec(:) INTEGER , INTENT(IN ) :: idatec(4) ! INTENT(IN ) :: idatec(:) INTEGER , INTENT(IN ) :: idatep(4) REAL(KIND=r8), INTENT(IN ) :: DELTAIN LOGICAL , INTENT(IN ) :: start INTEGER :: nCols INTEGER :: jb INTEGER :: ib INTEGER :: i INTEGER :: ios REAL(KIND=r8) :: rtime REAL(KIND=r8) :: time REAL(KIND=r8) :: time2 REAL(KIND=r8) :: jday REAL(KIND=r8) :: orbit REAL(KIND=r8) :: angle REAL(KIND=r8) :: xdecl REAL(KIND=r8) :: frac REAL(KIND=r8) :: sw REAL(KIND=r8) :: year REAL(KIND=r8) :: decmax REAL(KIND=r8) :: sols REAL(KIND=r8) :: season REAL(KIND=r8) :: day REAL(KIND=r8) :: dec REAL(KIND=r8) :: rfd REAL(KIND=r8) :: sind REAL(KIND=r8) :: cosd REAL(KIND=r8) :: hac REAL(KIND=r8) :: h REAL(KIND=r8) :: sr REAL(KIND=r8) :: ss REAL(KIND=r8) :: dawn REAL(KIND=r8) :: dusk REAL(KIND=r8) :: coshr REAL(KIND=r8) :: sunang REAL(KIND=r8) :: trans REAL(KIND=r8) :: btime REAL(KIND=r8) :: ctime REAL(KIND=r8) :: atime REAL(KIND=r8) :: fdiffuse REAL(KIND=r8) :: xlati (jbMax) REAL(KIND=r8) :: coszen(ibMax,jbMax) REAL(KIND=r8) :: cloud(ibMax,jbMax) REAL(KIND=r8) :: solad(ibMax,nband) REAL(KIND=r8) :: solai(ibMax,nband) REAL(KIND=r8) :: var(ibMax,jbMax) CHARACTER(LEN=255) :: namec CHARACTER(LEN=255) :: namep INTEGER :: open_status INTEGER :: reclen INTEGER :: irec INTEGER :: j INTEGER, PARAMETER :: input_unit=200 REAL(KIND=r8) , PARAMETER :: yrl=365.2500_r8 INTEGER, PARAMETER :: ok = 0 REAL(KIND=r8) , PARAMETER :: alon = 0.0_r8 REAL (KIND=r8), PARAMETER :: pie = 3.1415926e0_r8! constant pi=3.1415926e0 REAL (KIND=r8), PARAMETER :: pai12 = pie/12.0_r8 ! pi/12 LOGICAL :: here CHARACTER(LEN=255) :: PREFIX='GCM2IBIS' CHARACTER(LEN=255) :: SUFFIX='.dat' INTEGER :: iyeard INTEGER :: imonth INTEGER :: juliand REAL(KIND=r8) :: frachours REAL(KIND=r8) :: hours REAL(KIND=r8) :: xua REAL(KIND=r8) :: xta REAL(KIND=r8) :: xprec REAL(KIND=r8) :: xsin REAL(KIND=r8) :: xlin REAL(KIND=r8) :: xur REAL(KIND=r8) :: xes REAL(KIND=r8) :: xe REAL(KIND=r8) :: xqa REAL(KIND=r8) :: xpslc REAL(KIND=r8) :: xlh REAL(KIND=r8) :: xsh REAL(KIND=r8) :: sdelt REAL(KIND=r8) :: ratio REAL(KIND=r8) :: etime REAL(KIND=r8) :: xday REAL(KIND=r8) :: sindel REAL(KIND=r8) :: cosdel REAL(KIND=r8) :: fimxi REAL(KIND=r8) :: frh REAL(KIND=r8) :: cos2 (ibMax) REAL(KIND=r8) :: zenith1 (ibMax) REAL(KIND=r8) :: zenith2 (ibMax) REAL(KIND=r8) :: lonrad(ibMax) REAL(KIND=r8) :: swradtot (ibMax,jbMax) REAL(KIND=r8) :: diffswradtot(ibMax,jbMax) open_status=ok nCols=ibMax time=tod IF(MOD(tod,DELTAIN) == 0.0_r8 .OR. jdt == 1)THEN step=0 DO i=1,ibMax WRITE(namec,'(A12,i3.3)')'Clim_input.S',i INQUIRE (FILE='./inputdata/'//TRIM(namec),EXIST=here) IF (.NOT.here) THEN WRITE(0,*)'FILE ','./inputdata/'//TRIM(namec), 'Not Exist';STOP ELSE IF(OPENFILE)OPEN(input_unit+i,FILE='./inputdata/'//TRIM(namec),ACCESS='SEQUENTIAL',& FORM='FORMATTED',STATUS='OLD',ACTION='READ', & IOSTAT=open_status) IF( open_status == ok ) THEN IF(TRIM(icn_data) == 'ERA5')THEN READ(input_unit+i,'(3I8.8,9F8.3)') iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xlh,xsh ELSE IF(TRIM(icn_data) == 'GLDAS')THEN READ(input_unit+i,*) iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xpslc,xlh,xsh ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF DO j=1,jbMax latic(i,j) = lat_site(i)!degree longc(i,j) = lon_site(i)!degree t02mc(i,j) = xta+273.16_r8!K IF(TRIM(icn_data) == 'ERA5')THEN xpslc =( 1.15_r8*287.16*(xta+273.16_r8) )/100.0_r8 ! (rho*R*T) hPa xes=6.112*exp(17.67*xta/(243.5+xta))!hPa xe=(xur/100.0_r8)*xes !hPa xqa =xe*622.0/(xpslc-xe) !g/kg q02mc(i,j) = xqa/1000.0_r8!kg/kg ELSE IF(TRIM(icn_data) == 'GLDAS')THEN q02mc(i,j) = xur!kg/kg ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF uvesc(i,j) = xua!m/s vvesc(i,j) = xua!m/s pslcc(i,j) = xpslc!hPa IF(TRIM(icn_data) == 'ERA5')THEN precc(i,j) = xprec*86400/DELTAIN ocisc(i,j) = xsin!w/m2 olisc(i,j) = xlin!w/m2 slhfc(i,j) = xlh sshfc(i,j) = xsh ELSE IF(TRIM(icn_data) == 'GLDAS')THEN precc(i,j) = xprec ocisc(i,j) = xsin olisc(i,j) = xlin slhfc(i,j) = xlh sshfc(i,j) = xsh ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF END DO IF(.FALSE.)THEN WRITE(0,*)'latic',MAXVAL(latic),MINVAL(latic) WRITE(0,*)'longc',MAXVAL(longc),MINVAL(longc) WRITE(0,*)'t02mc',MAXVAL(t02mc),MINVAL(t02mc) WRITE(0,*)'q02mc',MAXVAL(q02mc),MINVAL(q02mc) WRITE(0,*)'uvesc',MAXVAL(uvesc),MINVAL(uvesc) WRITE(0,*)'vvesc',MAXVAL(vvesc),MINVAL(vvesc) WRITE(0,*)'pslcc',MAXVAL(pslcc),MINVAL(pslcc) WRITE(0,*)'precc',MAXVAL(precc),MINVAL(precc) WRITE(0,*)'ocisc',MAXVAL(ocisc),MINVAL(ocisc) WRITE(0,*)'olisc',MAXVAL(olisc),MINVAL(olisc) WRITE(0,*)'slhfc',MAXVAL(slhfc),MINVAL(slhfc) WRITE(0,*)'sshfc',MAXVAL(sshfc),MINVAL(sshfc) END IF ELSE WRITE(0,*)"Unable to OPEN file";STOP END IF END IF IF(start)REWIND(input_unit+i) END DO DO i=1,ibMax WRITE(namep,'(A13,i3.3)')'Clima_input.S',i INQUIRE (FILE='./inputdata/'//TRIM(namep),EXIST=here) IF (.NOT.here) THEN WRITE(0,*)'FILE ','./inputdata/'//TRIM(namep), 'Not Exist';STOP ELSE INQUIRE(IOLENGTH=reclen)latic IF(OPENFILE)OPEN(input_unit+ibMax+i,FILE='./inputdata/'//TRIM(namep),ACCESS='SEQUENTIAL',& FORM='FORMATTED',STATUS='OLD',ACTION='READ', & IOSTAT=open_status) IF( open_status == ok ) THEN IF(TRIM(icn_data) == 'ERA5')THEN READ(input_unit+ibMax+i,'(3I8.8,9F8.3)') iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xlh,xsh ELSE IF(TRIM(icn_data) == 'GLDAS')THEN IF(OPENFILE)READ(input_unit+ibMax+i,*) iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xpslc,xlh,xsh READ(input_unit+ibMax+i,*) iyeard,imonth,juliand,hours,xua,xta,xprec,xsin,xlin,xur,xpslc,xlh,xsh ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF DO j=1,jbMax latip(i,j) = lat_site(i)!degree longp(i,j) = lon_site(i)!degree t02mp(i,j) = xta+273.16_r8!K IF(TRIM(icn_data) == 'ERA5')THEN xpslc = (1.15_r8*287.16*(xta+273.16_r8))/100.0_r8 ! (rho*R*T) hPa xes=6.112*exp(17.67*xta/(243.5+xta))!hPa xe=(xur/100.0_r8)*xes !hPa xqa =xe*622.0/(xpslc-xe) !g/kg q02mp(i,j) = xqa/1000.0_r8!kg/kg ELSE IF(TRIM(icn_data) == 'GLDAS')THEN q02mc(i,j) = xur!kg/kg ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF uvesp(i,j) = xua!m/s vvesp(i,j) = xua!m/s pslcp(i,j) = xpslc!hPaxpslc!hPa IF(TRIM(icn_data) == 'ERA5')THEN precp(i,j) = xprec*86400/DELTAIN ocisp(i,j) = xsin!w/m2 olisp(i,j) = xlin!w/m2 slhfp(i,j) = xlh sshfp(i,j) = xsh ELSE IF(TRIM(icn_data) == 'GLDAS')THEN precp(i,j) = xprec ocisp(i,j) = xsin olisp(i,j) = xlin slhfp(i,j) = xlh sshfp(i,j) = xsh ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF END DO IF(.FALSE.)THEN WRITE(0,*)'latip',MAXVAL(latip),MINVAL(latip) WRITE(0,*)'longp',MAXVAL(longp),MINVAL(longp) WRITE(0,*)'t02mp',MAXVAL(t02mp),MINVAL(t02mp) WRITE(0,*)'q02mp',MAXVAL(q02mp),MINVAL(q02mp) WRITE(0,*)'uvesp',MAXVAL(uvesp),MINVAL(uvesp) WRITE(0,*)'vvesp',MAXVAL(vvesp),MINVAL(vvesp) WRITE(0,*)'pslcp',MAXVAL(pslcp),MINVAL(pslcp) WRITE(0,*)'precp',MAXVAL(precp),MINVAL(precp) WRITE(0,*)'ocisp',MAXVAL(ocisp),MINVAL(ocisp) WRITE(0,*)'olisp',MAXVAL(olisp),MINVAL(olisp) WRITE(0,*)'slhfp',MAXVAL(slhfp),MINVAL(slhfp) WRITE(0,*)'sshfp',MAXVAL(sshfp),MINVAL(sshfp) END IF ELSE WRITE(0,*)"Unable to OPEN file";STOP END IF END IF IF(start)REWIND(input_unit+ibMax+i) END DO WRITE(0,*)'READ ',TRIM(namec),' and ',TRIM(namep),iyeard,imonth END IF OPENFILE=.FALSE. DO jb=1,jbMax !pi --- 180 ! y x ! ! y = X*pi/180 !radianos ! lati (jbMax+1-jb) = (latic(1,jb)+90.0_r8)*pi/180.0_r8 xlati (jbMax+1-jb) = (latic(1,jb))*pi/180.0_r8 DO ib=1,nCols lonrad(ib) = longc(ib,1) gu (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*uvesc(ib,jb) & + ((step*dt)/(DELTAIN))*uvesp(ib,jb) gv (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*vvesc(ib,jb) & + ((step*dt)/(DELTAIN))*vvesp(ib,jb) gt (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) gq (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) & + ((step*dt)/(DELTAIN))*q02mp(ib,jb) gps (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*pslcc(ib,jb) & + ((step*dt)/(DELTAIN))*pslcp(ib,jb) var (ib, jb) = ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb) & + ((step*dt)/(DELTAIN))*precp(ib,jb) IF(TRIM(icn_data) == 'ERA5')THEN dcupr_sib (ib,jbMax+1-jb) = (var(ib,jb)/86400.0_r8)*1.0e+3_r8 ! PRINT*, var (ib, jb), ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb), ((step*dt)/(DELTAIN))*precp(ib,jb) dlspr_sib (ib,jbMax+1-jb) = 0.0_r8 ELSE IF(TRIM(icn_data) == 'GLDAS')THEN dcupr_sib (ib,jbMax+1-jb) = var(ib,jb)*1.0e+3_r8!/ convert m to mm ! PRINT*, var (ib, jb), ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb), ((step*dt)/(DELTAIN))*precp(ib,jb) dlspr_sib (ib,jbMax+1-jb) = 0.0_r8 ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF swdown (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*ocisc(ib,jb) & + ((step*dt)/(DELTAIN))*ocisp(ib,jb) dlwbot_sib(ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) & + ((step*dt)/(DELTAIN))*olisp(ib,jb) sigki = 1.0_r8 delsig = 0.010000_r8 tseam (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) tsea (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) slrad (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) & + ((step*dt)/(DELTAIN))*olisp(ib,jb) slhf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*slhfc(ib,jb) & + ((step*dt)/(DELTAIN))*slhfp(ib,jb) sshf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*sshf(ib,jb) & + ((step*dt)/(DELTAIN))*sshfp(ib,jb) ! zorl........zorl (i)= 100.0 *zgrav*speedm(i)*rhi(i) ! zgrav =0.032 /grav zorl (ib,jbMax+1-jb) = 100.0_r8*(0.32_r8/9.8_r8)*(sqrt( gu(ib,jbMax+1-jb)**2 + gv(ib,jbMax+1-jb)**2))*0.01 tsurf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) qsurf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) & + ((step*dt)/(DELTAIN))*q02mp(ib,jb) tmtx (ib,:,jbMax+1-jb) = 0.010000_r8 qmtx (ib,:,jbMax+1-jb) = 0.010000_r8 umtx (ib,:,jbMax+1-jb) = 0.010000_r8 IF(.TRUE.)THEN WRITE(0,*)'latip',MAXVAL(latip),MINVAL(latip) WRITE(0,*)'longp',MAXVAL(longp),MINVAL(longp) WRITE(0,*)'t02mp',MAXVAL(gt),MINVAL(gt) WRITE(0,*)'q02mp',MAXVAL(gq),MINVAL(gq) WRITE(0,*)'uvesp',MAXVAL(gu),MINVAL(gu) WRITE(0,*)'vvesp',MAXVAL(gv),MINVAL(gv) WRITE(0,*)'pslcp',MAXVAL(gps),MINVAL(gps) WRITE(0,*)'precp',MAXVAL(dcupr_sib),MINVAL(dcupr_sib) WRITE(0,*)'ocisp',MAXVAL(swdown),MINVAL(swdown) WRITE(0,*)'olisp',MAXVAL(dlwbot_sib),MINVAL(dlwbot_sib) WRITE(0,*)'slhfp',MAXVAL(slhfp),MINVAL(slhfp) WRITE(0,*)'sshfp',MAXVAL(sshfp),MINVAL(sshfp) END IF END DO END DO step=step+1 ! ! mon is the month used for vegetation data input ! DO jb=1,jbMax !DO ib=1,ibMax,1 ! month(ib,jb)=idatec(2) ! IF((((lati(jb)*180.0_r8)/3.1415926e0_r8)-90.0_r8) > 0.0_r8 ) THEN ! month(ib,jb) = month(ib,jb) + 6 ! IF(month(ib,jb).GE.13) month(ib,jb) = month(ib,jb)-12 ! END IF !END DO ! ! computation of astronomical parameters ! sdelt ;solar inclination ! etime ;correction factor to local time ! ratio ;factor relating to the distance between the earth and the sun ! CALL radtim(idatec,sdelt ,ratio ,etime ,tod ,xday ,yrl) sw=1370.0_r8*ratio sindel = SIN(sdelt) cosdel = COS(sdelt) fimxi = 24.0e0_r8 /360.0_r8 ctime = alon/15.0e0_r8 cos2 = 0.0e0_r8 frh = ( MOD(tod+0.03125_r8,3600.0_r8)-0.03125_r8)/3600.0_r8 DO ib=1,ibMax zenith1(ib) = sindel*COS(lati(jb)) ENDDO DO ib=1,ibMax btime = fimxi*lonrad(ib)+ctime atime = etime+pai12*(12.0_r8-idatec(1)-frh-btime) zenith2 (ib) = cosdel*SIN(lati(jb))*COS(atime) coszen (ib,jb) = MAX (0.0_r8,zenith1(ib) + zenith2(ib)) sunang = MAX( 0.01_r8, coszen (ib,jb) ) cloud(ib,jb) = (sw * sunang - swdown(ib,jb)) / (sw * sunang) cloud(ib,jb) = MAX(cloud(ib,jb),0.0_r8) cloud(ib,jb) = MIN(cloud(ib,jb),1.0_r8) !cloud(ib,jb) = MAX(0.58_r8 ,cloud(ib,jb)) END DO END DO ! ! do for each waveband ! DO ib = 1, nband DO jb=1,jbMax DO i=1,nCols ! calculate the solar transmission through the atmosphere ! using simple linear function of tranmission and cloud cover ! ! note that the 'cloud cover' data is typically obtained from ! sunshine hours -- not direct cloud observations ! ! where, cloud cover = 1 - sunshine fraction ! ! different authors present different values for the slope and ! intercept terms of this equation ! ! Friend, A: Parameterization of a global daily weather generator for ! terrestrial ecosystem and biogeochemical modelling, Ecological ! Modelling ! ! Spitters et al., 1986: Separating the diffuse and direct component ! of global radiation and its implications for modeling canopy ! photosynthesis, Part I: Components of incoming radiation, ! Agricultural and Forest Meteorology, 38, 217-229. ! ! A. Friend : trans = 0.251 + 0.509 * (1.0 - cloud(i)) ! Spitters et al. : trans = 0.200 + 0.560 * (1.0 - cloud(i)) ! ! we are using the values from A. Friend ! trans = 0.251_r8 + 0.509_r8 * (1.0_r8 - cloud(i,jb)) ! ! calculate the fraction of indirect (diffuse) solar radiation ! based upon the cloud cover ! ! note that these relationships typically are measured for either ! monthly or daily timescales, and may not be exactly appropriate ! for hourly calculations -- however, in ibis, cloud cover is fixed ! through the entire day so it may not make much difference ! ! method i -- ! ! we use a simple empirical relationships from Nikolov and Zeller (1992) ! ! Nikolov, N. and K.F. Zeller, 1992: A solar radiation algorithm for ecosystem ! dynamics models, Ecological Modelling, 61, 149-168. ! fdiffuse = 1.0045_r8 + 0.0435_r8 * trans - 3.5227_r8 * trans**2 + 2.6313_r8 * trans**3 ! IF (trans.GT.0.75_r8) fdiffuse = 0.166_r8 ! ! method ii -- ! ! another method was suggested by Spitters et al. (1986) based on ! long-term data from the Netherlands ! ! Spitters et al., 1986: Separating the diffuse and direct component ! of global radiation and its implications for modeling canopy ! photosynthesis, Part I: Components of incoming radiation, ! Agricultural and Forest Meteorology, 38, 217-229. ! ! if ((trans.eq.0.00).and.(trans.lt.0.07)) then ! fdiffuse = 1.0 ! else if ((trans.ge.0.07).and.(trans.lt.0.35)) then ! fdiffuse = 1.0 - 2.3 * (trans - 0.07)**2 ! else if ((trans.ge.0.35).and.(trans.lt.0.75)) then ! fdiffuse = 1.33 - 1.46 * trans ! else ! fdiffuse = 0.23 ! endif ! ! do for each waveband ! ! ! calculate the fraction in each waveband ! ! frac = 0.46 + 0.08 * float(ib - 1) frac = 0.427_r8 + 0.146_r8 * REAL((ib - 1),KIND=r8) ! ! calculate the direct and indirect solar radiation ! solad(i,ib) = swdown(i,jb) * frac * (1.0_r8 - fdiffuse) ! solai(i,ib) = swdown(i,jb) * frac * fdiffuse cosz_sib (i,jb) =coszen (i,jb) radvbc_sib (i,jb) =solad (i,1) radnbc_sib (i,jb) =solad (i,2) radvdc_sib (i,jb) =solai (i,1) radndc_sib (i,jb) =solai (i,2) ! END DO END DO END DO DO jb=1,jbMax DO i=1,nCols swradtot(i,jb) =radvbc_sib (i,jb)+ radnbc_sib (i,jb)+ radvdc_sib (i,jb)+ radndc_sib (i,jb) diffswradtot(i,jb)=swdown(i,jb)-swradtot(i,jb) END DO END DO DO jb=1,jbMax DO i=1,nCols radvbc_sib (i,jb) =radvbc_sib (i,jb) + 0.25_r8*diffswradtot(i,jb) radnbc_sib (i,jb) =radnbc_sib (i,jb) + 0.25_r8*diffswradtot(i,jb) radvdc_sib (i,jb) =radvdc_sib (i,jb) + 0.25_r8*diffswradtot(i,jb) radndc_sib (i,jb) =radndc_sib (i,jb) + 0.25_r8*diffswradtot(i,jb) END DO END DO RETURN 100 STOP END SUBROUTINE RunAtmosModel1D SUBROUTINE RunAtmosModel2D(rec,jdt,tod,nband,dt,idate,idatec,idatep,DELTAIN,start,iMask) IMPLICIT NONE INTEGER , INTENT(INOUT) :: rec REAL(KIND=r8), INTENT(IN ) :: tod INTEGER , INTENT(IN ) :: jdt INTEGER , INTENT(IN ) :: nband ! INTENT(IN ) :: nband REAL(KIND=r8), INTENT(IN ) :: dt ! INTENT(IN ) :: dt INTEGER , INTENT(IN ) :: idate (4) ! INTENT(IN ) :: idatec(:) INTEGER , INTENT(IN ) :: idatec(4) ! INTENT(IN ) :: idatec(:) INTEGER , INTENT(IN ) :: idatep(4) REAL(KIND=r8), INTENT(IN ) :: DELTAIN LOGICAL , INTENT(IN ) :: start INTEGER(KIND=i8) , INTENT(IN ) :: imask (ibMax, jbMax) INTEGER :: nCols INTEGER :: jb INTEGER :: ib INTEGER :: i INTEGER :: ios REAL(KIND=r8) :: rtime REAL(KIND=r8) :: time REAL(KIND=r8) :: time2 REAL(KIND=r8) :: jday REAL(KIND=r8) :: orbit REAL(KIND=r8) :: angle REAL(KIND=r8) :: xdecl REAL(KIND=r8) :: frac REAL(KIND=r8) :: sw REAL(KIND=r8) :: year REAL(KIND=r8) :: decmax REAL(KIND=r8) :: sols REAL(KIND=r8) :: season REAL(KIND=r8) :: day REAL(KIND=r8) :: dec REAL(KIND=r8) :: rfd REAL(KIND=r8) :: sind REAL(KIND=r8) :: cosd REAL(KIND=r8) :: hac REAL(KIND=r8) :: h REAL(KIND=r8) :: sr REAL(KIND=r8) :: ss REAL(KIND=r8) :: dawn REAL(KIND=r8) :: dusk REAL(KIND=r8) :: coshr REAL(KIND=r8) :: sunang REAL(KIND=r8) :: trans REAL(KIND=r8) :: btime REAL(KIND=r8) :: ctime REAL(KIND=r8) :: atime REAL(KIND=r8) :: fdiffuse REAL(KIND=r8) :: xlati (jbMax) REAL(KIND=r8) :: coszen(ibMax,jbMax) REAL(KIND=r8) :: cloud(ibMax,jbMax) REAL(KIND=r8) :: solad(ibMax,nband) REAL(KIND=r8) :: solai(ibMax,nband) REAL(KIND=r8) :: var(ibMax,jbMax) CHARACTER(LEN=255) :: namec CHARACTER(LEN=255) :: namep INTEGER :: open_status INTEGER :: reclen INTEGER :: irec INTEGER :: j INTEGER, PARAMETER :: input_unit=200 REAL(KIND=r8) , PARAMETER :: yrl=365.2500_r8 INTEGER, PARAMETER :: ok = 0 REAL(KIND=r8) , PARAMETER :: alon = 0.0_r8 REAL (KIND=r8), PARAMETER :: pie = 3.1415926e0_r8! constant pi=3.1415926e0 REAL (KIND=r8), PARAMETER :: pai12 = pie/12.0_r8 ! pi/12 LOGICAL :: here CHARACTER(LEN=255) :: PREFIX='GCM2IBIS' CHARACTER(LEN=255) :: SUFFIX='.dat' INTEGER :: iyeard INTEGER :: imonth INTEGER :: juliand REAL(KIND=r8) :: frachours REAL(KIND=r8) :: hours REAL(KIND=r8) :: xua REAL(KIND=r8) :: xta REAL(KIND=r8) :: xprec REAL(KIND=r8) :: xsin REAL(KIND=r8) :: xlin REAL(KIND=r8) :: xur REAL(KIND=r8) :: xes REAL(KIND=r8) :: xe REAL(KIND=r8) :: xqa REAL(KIND=r8) :: xpslc REAL(KIND=r8) :: xlh=0 REAL(KIND=r8) :: xsh=0 REAL(KIND=r8) :: sdelt REAL(KIND=r8) :: ratio REAL(KIND=r8) :: etime REAL(KIND=r8) :: xday REAL(KIND=r8) :: sindel REAL(KIND=r8) :: cosdel REAL(KIND=r8) :: fimxi REAL(KIND=r8) :: frh REAL(KIND=r8) :: cos2 (ibMax) REAL(KIND=r8) :: zenith1 (ibMax) REAL(KIND=r8) :: zenith2 (ibMax) REAL(KIND=r8) :: lonrad(ibMax) REAL(KIND=r8) :: xua_stepc(ibMax,jbMax) REAL(KIND=r8) :: xta_stepc(ibMax,jbMax) REAL(KIND=r8) :: xpres_stepc(ibMax,jbMax) REAL(KIND=r8) :: xprec_stepc(ibMax,jbMax) REAL(KIND=r8) :: xsin_stepc(ibMax,jbMax) REAL(KIND=r8) :: xlin_stepc(ibMax,jbMax) REAL(KIND=r8) :: xur_stepc(ibMax,jbMax) REAL(KIND=r8) :: xqa_stepc(ibMax,jbMax) REAL(KIND=r8) :: xua_stepp(ibMax,jbMax) REAL(KIND=r8) :: xta_stepp(ibMax,jbMax) REAL(KIND=r8) :: xpres_stepp(ibMax,jbMax) REAL(KIND=r8) :: xprec_stepp(ibMax,jbMax) REAL(KIND=r8) :: xsin_stepp(ibMax,jbMax) REAL(KIND=r8) :: xlin_stepp(ibMax,jbMax) REAL(KIND=r8) :: xur_stepp(ibMax,jbMax) REAL(KIND=r8) :: xqa_stepp(ibMax,jbMax) REAL(KIND=r8) :: swradtot (ibMax,jbMax) REAL(KIND=r8) :: diffswradtot(ibMax,jbMax) open_status=ok nCols=ibMax time=tod IF(MOD(tod,DELTAIN) == 0.0_r8 .OR. jdt == 1)THEN step=0 rec=rec+1 IF(TRIM(icn_data) == 'ERA5')THEN !PK CALL RunReadNETCDF(rec,ncid_prec,latic,latip,longc,longp,xprec_stepc,xprec_stepp,nFile=1) !PK CALL RunReadNETCDF(rec,ncid_temp,latic,latip,longc,longp,xta_stepc ,xta_stepp ,nFile=2) !PK CALL RunReadNETCDF(rec,ncid_rhum,latic,latip,longc,longp,xur_stepc ,xur_stepp ,nFile=3) !PK CALL RunReadNETCDF(rec,ncid_wind,latic,latip,longc,longp,xua_stepc ,xua_stepp ,nFile=4) !PK CALL RunReadNETCDF(rec,ncid_swrd,latic,latip,longc,longp,xsin_stepc ,xsin_stepp ,nFile=5) !PK CALL RunReadNETCDF(rec,ncid_lwrd,latic,latip,longc,longp,xlin_stepc ,xlin_stepp, nFile=6) ELSE IF(TRIM(icn_data) == 'GLDAS')THEN !PK CALL RunReadNETCDF(rec,ncid_pres,latic,latip,longc,longp,xpres_stepc,xpres_stepp,nFile=1) !PK CALL RunReadNETCDF(rec,ncid_temp,latic,latip,longc,longp,xta_stepc ,xta_stepp ,nFile=2) !PK CALL RunReadNETCDF(rec,ncid_umes,latic,latip,longc,longp,xqa_stepc ,xqa_stepp ,nFile=3) !PK CALL RunReadNETCDF(rec,ncid_wind,latic,latip,longc,longp,xua_stepc ,xua_stepp ,nFile=4) !PK CALL RunReadNETCDF(rec,ncid_prec,latic,latip,longc,longp,xprec_stepc,xprec_stepp,nFile=5) !PK CALL RunReadNETCDF(rec,ncid_swrd,latic,latip,longc,longp,xsin_stepc ,xsin_stepp ,nFile=6) !PK CALL RunReadNETCDF(rec,ncid_lwrd,latic,latip,longc,longp,xlin_stepc ,xlin_stepp, nFile=7) ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF DO j=1,jbMax DO i=1,ibMax !IF (iMask(i,jbMax+1-j) >= 1) THEN t02mc(i,j) = xta_stepc(i,j)+273.16_r8 !K !xpslc = 997.0_r8 !hPa IF(TRIM(icn_data) == 'ERA5')THEN xpslc = (1.15_r8*287.16*(xta_stepc(i,j)+273.16_r8) )/100.0_r8 ! (rho*R*T) hPa xes=6.112*exp(17.67*xta_stepc(i,j)/(243.5+xta_stepc(i,j)))!hPa xe=(xur_stepc(i,j)/100.0_r8)*xes !hPa xqa =xe*622.0/(xpslc-xe) !g/kg q02mc(i,j) = xqa/1000.0_r8 !kg/kg ELSE IF(TRIM(icn_data) == 'GLDAS')THEN xpslc = (1.15_r8*287.16*(xta_stepc(i,j)+273.16_r8) )/100.0_r8 ! (rho*R*T) hPa q02mc(i,j) = xqa_stepc(i,j) ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF uvesc(i,j) = xua_stepc(i,j) !m/s vvesc(i,j) = xua_stepc(i,j) !m/s IF(TRIM(icn_data) == 'ERA5')THEN pslcc(i,j) = xpslc !hPa ELSE IF(TRIM(icn_data) == 'GLDAS')THEN IF(xpres_stepc(i,j) < 100) THEN pslcc(i,j) =xpslc ELSE pslcc(i,j) = xpres_stepc(i,j) END IF ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF IF(TRIM(icn_data) == 'ERA5')THEN precc(i,j) = xprec_stepc(i,j)*86400/DELTAIN !mm/dia ELSE IF(TRIM(icn_data) == 'GLDAS')THEN precc(i,j) = xprec_stepc(i,j) !*86400/DELTAIN !mm/dia ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF ocisc(i,j) = xsin_stepc(i,j) !w/m2 olisc(i,j) = xlin_stepc(i,j) !w/m2 ! END IF END DO END DO IF(.FALSE.)THEN WRITE(0,*)'latic',MAXVAL(latic),MINVAL(latic) WRITE(0,*)'longc',MAXVAL(longc),MINVAL(longc) WRITE(0,*)'t02mc',MAXVAL(t02mc),MINVAL(t02mc) WRITE(0,*)'q02mc',MAXVAL(q02mc),MINVAL(q02mc) WRITE(0,*)'uvesc',MAXVAL(uvesc),MINVAL(uvesc) WRITE(0,*)'vvesc',MAXVAL(vvesc),MINVAL(vvesc) WRITE(0,*)'pslcc',MAXVAL(pslcc),MINVAL(pslcc) WRITE(0,*)'precc',MAXVAL(precc),MINVAL(precc) WRITE(0,*)'ocisc',MAXVAL(ocisc),MINVAL(ocisc) WRITE(0,*)'olisc',MAXVAL(olisc),MINVAL(olisc) WRITE(0,*)'slhfc',MAXVAL(slhfc),MINVAL(slhfc) WRITE(0,*)'sshfc',MAXVAL(sshfc),MINVAL(sshfc) END IF DO j=1,jbMax DO i=1,ibMax !IF (iMask(i,jbMax+1-j) >= 1) THEN t02mp(i,j) = xta_stepp(i,j)+273.16_r8 !K !xpslc = 1000.0_r8 !hPa IF(TRIM(icn_data) == 'ERA5')THEN xpslc = (1.15_r8*287.16*(xta_stepp(i,j)+273.16_r8) )/100.0_r8 ! (rho*R*T) hPa xes = 6.112*exp(17.67*xta_stepp(i,j)/(243.5+xta_stepp(i,j)))!hPa xe = (xur_stepp(i,j)/100.0_r8)*xes !hPa xqa = xe*622.0/(xpslc-xe) !g/kg q02mp(i,j) = xqa/1000.0_r8 !kg/kg ELSE IF(TRIM(icn_data) == 'GLDAS')THEN xpslc = (1.15_r8*287.16*(xta_stepp(i,j)+273.16_r8) )/100.0_r8 ! (rho*R*T) hPa q02mp(i,j) = xqa_stepp(i,j) ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF uvesp(i,j) = xua_stepp(i,j) !m/s vvesp(i,j) = xua_stepp(i,j) !m/s IF(TRIM(icn_data) == 'ERA5')THEN pslcp(i,j) = xpslc !hPa ELSE IF(TRIM(icn_data) == 'GLDAS')THEN IF(xpres_stepp(i,j) < 100) THEN pslcp(i,j) =xpslc ELSE pslcp(i,j) = xpres_stepp(i,j) END IF ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF IF(TRIM(icn_data) == 'ERA5')THEN precp(i,j) = xprec_stepp(i,j)*86400/DELTAIN !mm/dia ELSE IF(TRIM(icn_data) == 'GLDAS')THEN precp(i,j) = xprec_stepp(i,j)!*86400/DELTAIN !mm/dia ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF ocisp(i,j) = xsin_stepp(i,j) !w/m2 olisp(i,j) = xlin_stepp(i,j) !w/m2 slhfp(i,j) = xlh sshfp(i,j) = xsh ! END IF END DO END DO IF(.FALSE.)THEN WRITE(0,*)'latip',MAXVAL(latip),MINVAL(latip) WRITE(0,*)'longp',MAXVAL(longp),MINVAL(longp) WRITE(0,*)'t02mp',MAXVAL(t02mp),MINVAL(t02mp) WRITE(0,*)'q02mp',MAXVAL(q02mp),MINVAL(q02mp) WRITE(0,*)'uvesp',MAXVAL(uvesp),MINVAL(uvesp) WRITE(0,*)'vvesp',MAXVAL(vvesp),MINVAL(vvesp) WRITE(0,*)'pslcp',MAXVAL(pslcp),MINVAL(pslcp) WRITE(0,*)'precp',MAXVAL(precp),MINVAL(precp) WRITE(0,*)'ocisp',MAXVAL(ocisp),MINVAL(ocisp) WRITE(0,*)'olisp',MAXVAL(olisp),MINVAL(olisp) WRITE(0,*)'slhfp',MAXVAL(slhfp),MINVAL(slhfp) WRITE(0,*)'sshfp',MAXVAL(sshfp),MINVAL(sshfp) END IF END IF DO jb=1,jbMax DO ib=1,nCols ! IF (iMask(ib,jb) >= 1) THEN gu (ib,jbMax+1-jb) = 0.00000_r8 gv (ib,jbMax+1-jb) = 0.00000_r8 gt (ib,jbMax+1-jb) = 0.00000_r8 gq (ib,jbMax+1-jb) = 0.00000_r8 gps (ib,jbMax+1-jb) = 0.00000_r8 var (ib, jb) = 0.00000_r8 dcupr_sib (ib,jbMax+1-jb) = 0.00000_r8 dlspr_sib (ib,jbMax+1-jb) = 0.00000_r8 swdown (ib,jbMax+1-jb) = 0.00000_r8 dlwbot_sib(ib,jbMax+1-jb) = 0.00000_r8 sigki = 0.00000_r8 delsig = 0.00000_r8 tseam (ib,jbMax+1-jb) = 0.00000_r8 tsea (ib,jbMax+1-jb) = 0.00000_r8 slrad (ib,jbMax+1-jb) = 0.00000_r8 slhf (ib,jbMax+1-jb) = 0.00000_r8 sshf (ib,jbMax+1-jb) = 0.00000_r8 zorl (ib,jbMax+1-jb) = 0.00000_r8 tsurf (ib,jbMax+1-jb) = 0.00000_r8 qsurf (ib,jbMax+1-jb) = 0.00000_r8 tmtx (ib,:,jbMax+1-jb) = 0.00000_r8 qmtx (ib,:,jbMax+1-jb) = 0.00000_r8 umtx (ib,:,jbMax+1-jb) = 0.00000_r8 ! END IF END DO !pi --- 180 ! y x ! ! y = X*pi/180 !radianos ! lati (jbMax+1-jb) = (latic(1,jb)+90.0_r8)*pi/180.0_r8 xlati (jb) = (latic(1,jb)+90.0_r8)*pi/180.0_r8 DO ib=1,nCols lonrad(ib) = longc(ib,1) IF (iMask(ib,jbMax+1-jb) >= 1) THEN gu (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*uvesc(ib,jb) & + ((step*dt)/(DELTAIN))*uvesp(ib,jb) gv (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*vvesc(ib,jb) & + ((step*dt)/(DELTAIN))*vvesp(ib,jb) gt (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) gq (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) & + ((step*dt)/(DELTAIN))*q02mp(ib,jb) gps (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*pslcc(ib,jb) & + ((step*dt)/(DELTAIN))*pslcp(ib,jb) var (ib, jb) = ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb) & + ((step*dt)/(DELTAIN))*precp(ib,jb) IF(TRIM(icn_data) == 'ERA5')THEN dcupr_sib (ib,jbMax+1-jb) = (var(ib,jb)/86400.0_r8)*1.0e+3_r8 ! PRINT*, var (ib, jb), ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb), ((step*dt)/(DELTAIN))*precp(ib,jb) dlspr_sib (ib,jbMax+1-jb) = 0.0_r8 ELSE IF(TRIM(icn_data) == 'GLDAS')THEN dcupr_sib (ib,jbMax+1-jb) = var(ib,jb)*1.0e+3_r8!/86400.0_r8)*1.0e+3_r8 ! PRINT*, var (ib, jb), ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb), ((step*dt)/(DELTAIN))*precp(ib,jb) dlspr_sib (ib,jbMax+1-jb) = 0.0_r8 ELSE PRINT*,'ERROR',TRIM(icn_data),'diff','ERA5','GLDAS' STOP END IF swdown (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*ocisc(ib,jb) & + ((step*dt)/(DELTAIN))*ocisp(ib,jb) dlwbot_sib(ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) & + ((step*dt)/(DELTAIN))*olisp(ib,jb) sigki = 1.0_r8 delsig = 0.010000_r8 tseam (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) tsea (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) slrad (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) & + ((step*dt)/(DELTAIN))*olisp(ib,jb) slhf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*slhfc(ib,jb) & + ((step*dt)/(DELTAIN))*slhfp(ib,jb) sshf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*sshf(ib,jb) & + ((step*dt)/(DELTAIN))*sshfp(ib,jb) ! zorl........zorl (i)= 100.0 *zgrav*speedm(i)*rhi(i) ! zgrav =0.032 /grav zorl (ib,jbMax+1-jb) = 100.0_r8*(0.32_r8/9.8_r8)*(sqrt( gu(ib,jbMax+1-jb)**2 + gv(ib,jbMax+1-jb)**2))*0.01 tsurf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) qsurf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) & + ((step*dt)/(DELTAIN))*q02mp(ib,jb) tmtx (ib,:,jbMax+1-jb) = 0.010000_r8 qmtx (ib,:,jbMax+1-jb) = 0.010000_r8 umtx (ib,:,jbMax+1-jb) = 0.010000_r8 endif END DO END DO step=step+1 IF(.TRUE.)THEN WRITE(0,*)'latip',MAXVAL(latip),MINVAL(latip) WRITE(0,*)'longp',MAXVAL(longp),MINVAL(longp) WRITE(0,*)'t02mp',MAXVAL(gt),MINVAL(gt) WRITE(0,*)'q02mp',MAXVAL(gq),MINVAL(gq) WRITE(0,*)'uvesp',MAXVAL(gu),MINVAL(gu) WRITE(0,*)'vvesp',MAXVAL(gv),MINVAL(gv) WRITE(0,*)'pslcp',MAXVAL(gps),MINVAL(gps) WRITE(0,*)'precp',MAXVAL(dcupr_sib),MINVAL(dcupr_sib) WRITE(0,*)'ocisp',MAXVAL(swdown),MINVAL(swdown) WRITE(0,*)'olisp',MAXVAL(dlwbot_sib),MINVAL(dlwbot_sib) WRITE(0,*)'slhfp',MAXVAL(slhfp),MINVAL(slhfp) WRITE(0,*)'sshfp',MAXVAL(sshfp),MINVAL(sshfp) END IF ! ! mon is the month used for vegetation data input ! DO jb=1,jbMax !DO ib=1,ibMax,1 ! month(ib,jb)=idatec(2) ! IF((((lati(jb)*180.0_r8)/3.1415926e0_r8)-90.0_r8) > 0.0_r8 ) THEN ! month(ib,jb) = month(ib,jb) + 6 ! IF(month(ib,jb).GE.13) month(ib,jb) = month(ib,jb)-12 ! END IF !END DO ! ! computation of astronomical parameters ! sdelt ;solar inclination ! etime ;correction factor to local time ! ratio ;factor relating to the distance between the earth and the sun ! CALL radtim(idatec,sdelt ,ratio ,etime ,tod ,xday ,yrl) sw=1370.0_r8*ratio sindel = SIN(sdelt) cosdel = COS(sdelt) fimxi = 24.0e0_r8 /360.0_r8 ctime = alon/15.0e0_r8 cos2 = 0.0e0_r8 frh = ( MOD(tod+0.03125_r8,3600.0_r8)-0.03125_r8)/3600.0_r8 DO ib=1,ibMax zenith1(ib) = sindel*COS(xlati(jb)) ENDDO DO ib=1,ibMax btime = fimxi*lonrad(ib)+ctime atime = etime+pai12*(12.0_r8-idatec(1)-frh-btime) zenith2 (ib) = cosdel*SIN(xlati(jb))*COS(atime) coszen (ib,jb) = MAX (0.0_r8,zenith1(ib) + zenith2(ib)) sunang = MAX( 0.01_r8, coszen (ib,jb) ) cloud(ib,jb) = (sw * sunang - swdown(ib,jb)) / (sw * sunang) cloud(ib,jb) = MAX(cloud(ib,jb),0.0_r8) cloud(ib,jb) = MIN(cloud(ib,jb),1.0_r8) !cloud(ib,jb) = MAX(0.58_r8 ,cloud(ib,jb)) END DO END DO ! ! do for each waveband ! DO ib = 1, nband DO jb=1,jbMax DO i=1,nCols ! calculate the solar transmission through the atmosphere ! using simple linear function of tranmission and cloud cover ! ! note that the 'cloud cover' data is typically obtained from ! sunshine hours -- not direct cloud observations ! ! where, cloud cover = 1 - sunshine fraction ! ! different authors present different values for the slope and ! intercept terms of this equation ! ! Friend, A: Parameterization of a global daily weather generator for ! terrestrial ecosystem and biogeochemical modelling, Ecological ! Modelling ! ! Spitters et al., 1986: Separating the diffuse and direct component ! of global radiation and its implications for modeling canopy ! photosynthesis, Part I: Components of incoming radiation, ! Agricultural and Forest Meteorology, 38, 217-229. ! ! A. Friend : trans = 0.251 + 0.509 * (1.0 - cloud(i)) ! Spitters et al. : trans = 0.200 + 0.560 * (1.0 - cloud(i)) ! ! we are using the values from A. Friend ! trans = 0.251_r8 + 0.509_r8 * (1.0_r8 - cloud(i,jb)) ! ! calculate the fraction of indirect (diffuse) solar radiation ! based upon the cloud cover ! ! note that these relationships typically are measured for either ! monthly or daily timescales, and may not be exactly appropriate ! for hourly calculations -- however, in ibis, cloud cover is fixed ! through the entire day so it may not make much difference ! ! method i -- ! ! we use a simple empirical relationships from Nikolov and Zeller (1992) ! ! Nikolov, N. and K.F. Zeller, 1992: A solar radiation algorithm for ecosystem ! dynamics models, Ecological Modelling, 61, 149-168. ! fdiffuse = 1.0045_r8 + 0.0435_r8 * trans - 3.5227_r8 * trans**2 + 2.6313_r8 * trans**3 ! IF (trans.GT.0.75_r8) fdiffuse = 0.166_r8 ! ! method ii -- ! ! another method was suggested by Spitters et al. (1986) based on ! long-term data from the Netherlands ! ! Spitters et al., 1986: Separating the diffuse and direct component ! of global radiation and its implications for modeling canopy ! photosynthesis, Part I: Components of incoming radiation, ! Agricultural and Forest Meteorology, 38, 217-229. ! ! if ((trans.eq.0.00).and.(trans.lt.0.07)) then ! fdiffuse = 1.0 ! else if ((trans.ge.0.07).and.(trans.lt.0.35)) then ! fdiffuse = 1.0 - 2.3 * (trans - 0.07)**2 ! else if ((trans.ge.0.35).and.(trans.lt.0.75)) then ! fdiffuse = 1.33 - 1.46 * trans ! else ! fdiffuse = 0.23 ! endif ! ! do for each waveband ! ! ! calculate the fraction in each waveband ! ! frac = 0.46 + 0.08 * float(ib - 1) frac = 0.427_r8 + 0.146_r8 * REAL((ib - 1),KIND=r8) ! ! calculate the direct and indirect solar radiation ! solad(i,ib) = swdown(i,jb) * frac * (1.0_r8 - fdiffuse) ! solai(i,ib) = swdown(i,jb) * frac * fdiffuse cosz_sib (i,jb) =coszen (i,jb) radvbc_sib (i,jb) =solad (i,1) radnbc_sib (i,jb) =solad (i,2) radvdc_sib (i,jb) =solai (i,1) radndc_sib (i,jb) =solai (i,2) ! END DO END DO END DO DO jb=1,jbMax DO i=1,nCols swradtot(i,jb) =radvbc_sib (i,jb)+ radnbc_sib (i,jb)+ radvdc_sib (i,jb)+ radndc_sib (i,jb) diffswradtot(i,jb)=swdown(i,jb)-swradtot(i,jb) END DO END DO DO jb=1,jbMax DO i=1,nCols radvbc_sib (i,jb) =radvbc_sib (i,jb) + 0.25_r8*diffswradtot(i,jb) radnbc_sib (i,jb) =radnbc_sib (i,jb) + 0.25_r8*diffswradtot(i,jb) radvdc_sib (i,jb) =radvdc_sib (i,jb) + 0.25_r8*diffswradtot(i,jb) radndc_sib (i,jb) =radndc_sib (i,jb) + 0.25_r8*diffswradtot(i,jb) END DO END DO RETURN 100 STOP END SUBROUTINE RunAtmosModel2D SUBROUTINE RunAtmosModel(jdt,tod,nband,dt,idate,idatec,idatep,DELTAIN) IMPLICIT NONE REAL(KIND=r8), INTENT(IN ) :: tod INTEGER , INTENT(IN ) :: jdt INTEGER , INTENT(IN ) :: nband ! INTENT(IN ) :: nband REAL(KIND=r8), INTENT(IN ) :: dt ! INTENT(IN ) :: dt INTEGER , INTENT(IN ) :: idate (4) ! INTENT(IN ) :: idatec(:) INTEGER , INTENT(IN ) :: idatec(4) ! INTENT(IN ) :: idatec(:) INTEGER , INTENT(IN ) :: idatep(4) REAL(KIND=r8), INTENT(IN ) :: DELTAIN INTEGER :: nCols INTEGER :: jb INTEGER :: ib INTEGER :: i INTEGER :: ios REAL(KIND=r8) :: rtime REAL(KIND=r8) :: time REAL(KIND=r8) :: time2 REAL(KIND=r8) :: jday REAL(KIND=r8) :: orbit REAL(KIND=r8) :: angle REAL(KIND=r8) :: xdecl REAL(KIND=r8) :: frac REAL(KIND=r8) :: sw REAL(KIND=r8) :: year REAL(KIND=r8) :: decmax REAL(KIND=r8) :: sols REAL(KIND=r8) :: season REAL(KIND=r8) :: day REAL(KIND=r8) :: dec REAL(KIND=r8) :: rfd REAL(KIND=r8) :: sind REAL(KIND=r8) :: cosd REAL(KIND=r8) :: hac REAL(KIND=r8) :: h REAL(KIND=r8) :: sr REAL(KIND=r8) :: ss REAL(KIND=r8) :: dawn REAL(KIND=r8) :: dusk REAL(KIND=r8) :: coshr REAL(KIND=r8) :: sunang REAL(KIND=r8) :: trans REAL(KIND=r8) :: btime REAL(KIND=r8) :: ctime REAL(KIND=r8) :: atime REAL(KIND=r8) :: fdiffuse REAL(KIND=r8) :: xlati (jbMax) REAL(KIND=r8) :: coszen(ibMax,jbMax) REAL(KIND=r8) :: cloud(ibMax,jbMax) REAL(KIND=r8) :: swdown(ibMax,jbMax) REAL(KIND=r8) :: solad(ibMax,nband) REAL(KIND=r8) :: solai(ibMax,nband) REAL(KIND=r8) :: var(ibMax,jbMax) CHARACTER(LEN=255) :: namec CHARACTER(LEN=255) :: namep INTEGER :: open_status INTEGER :: reclen INTEGER :: irec INTEGER, PARAMETER :: input_unit=200 REAL(KIND=r8) , PARAMETER :: yrl=365.2500_r8 INTEGER, PARAMETER :: ok = 0 REAL(KIND=r8) , PARAMETER :: alon = 0.0_r8 REAL (KIND=r8), PARAMETER :: pie = 3.1415926e0_r8! constant pi=3.1415926e0 REAL (KIND=r8), PARAMETER :: pai12 = pie/12.0_r8 ! pi/12 LOGICAL :: here CHARACTER(LEN=255) :: PREFIX='GCM2IBIS' CHARACTER(LEN=255) :: SUFFIX='.dat' REAL(KIND=r8) :: sdelt REAL(KIND=r8) :: ratio REAL(KIND=r8) :: etime REAL(KIND=r8) :: xday REAL(KIND=r8) :: sindel REAL(KIND=r8) :: cosdel REAL(KIND=r8) :: fimxi REAL(KIND=r8) :: frh REAL(KIND=r8) :: cos2 (ibMax) REAL(KIND=r8) :: zenith1 (ibMax) REAL(KIND=r8) :: zenith2 (ibMax) REAL(KIND=r8) :: lonrad(ibMax) nCols=ibMax time=tod IF(MOD(tod,DELTAIN) == 0.0_r8 .OR. jdt == 1)THEN step=0 WRITE(namec,'(A8,2(I4.4,3i2.2),A4)')TRIM(PREFIX),& idate(4),idate(2),idate(3),idate(1),& idatec(4),idatec(2),idatec(3),idatec(1),& TRIM(SUFFIX) INQUIRE (FILE='./inputdata/'//TRIM(namec),EXIST=here) IF (.NOT.here) THEN WRITE(0,*)'FILE ','./inputdata/'//TRIM(namec), 'Not Exist';STOP ELSE INQUIRE(IOLENGTH=reclen)latic OPEN(input_unit,FILE='./inputdata/'//TRIM(namec),ACCESS='DIRECT',& FORM='UNFORMATTED',STATUS='OLD',RECL=reclen, & IOSTAT=open_status) IF( open_status == ok ) THEN irec=1 READ(input_unit,rec=irec,iostat = ios)latic;IF(ios /= 0) STOP WRITE(0,*)'latic',MAXVAL(latic),MINVAL(latic) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)longc;IF(ios /= 0) STOP WRITE(0,*)'longc',MAXVAL(longc),MINVAL(longc) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)t02mc;IF(ios /= 0) STOP WRITE(0,*)'t02mc',MAXVAL(t02mc),MINVAL(t02mc) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)q02mc;IF(ios /= 0) STOP WRITE(0,*)'q02mc',MAXVAL(q02mc),MINVAL(q02mc) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)uvesc;IF(ios /= 0) STOP WRITE(0,*)'uvesc',MAXVAL(uvesc),MINVAL(uvesc) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)vvesc;IF(ios /= 0) STOP WRITE(0,*)'vvesc',MAXVAL(vvesc),MINVAL(vvesc) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)pslcc;IF(ios /= 0) STOP WRITE(0,*)'pslcc',MAXVAL(pslcc),MINVAL(pslcc) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)precc;IF(ios /= 0) STOP WRITE(0,*)'precc',MAXVAL(precc),MINVAL(precc) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)ocisc;IF(ios /= 0) STOP WRITE(0,*)'ocisc',MAXVAL(ocisc),MINVAL(ocisc) irec=irec+1 READ(input_unit,rec=irec,iostat = ios)olisc;IF(ios /= 0) STOP WRITE(0,*)'olisc',MAXVAL(olisc),MINVAL(olisc) CLOSE(input_unit,STATUS='KEEP') ELSE WRITE(0,*)"Unable to OPEN file";STOP END IF END IF WRITE(namep,'(A8,2(I4.4,3i2.2),A4)')TRIM(PREFIX),& idate(4),idate(2),idate(3),idate(1),& idatep(4),idatep(2),idatep(3),idatep(1),& TRIM(SUFFIX) INQUIRE (FILE='./inputdata/'//TRIM(namep),EXIST=here) IF (.NOT.here) THEN WRITE(0,*)'FILE ','./inputdata/'//TRIM(namep), 'Not Exist';STOP ELSE INQUIRE(IOLENGTH=reclen)latic OPEN(input_unit,FILE='./inputdata/'//TRIM(namep),ACCESS='DIRECT',& FORM='UNFORMATTED',STATUS='OLD',RECL=reclen, & IOSTAT=open_status) IF( open_status == ok ) THEN irec=1 READ(input_unit,rec=irec,iostat = ios)latip;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)longp;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)t02mp;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)q02mp;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)uvesp;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)vvesp;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)pslcp;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)precp;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)ocisp;IF(ios /= 0) STOP irec=irec+1 READ(input_unit,rec=irec,iostat = ios)olisp;IF(ios /= 0) STOP CLOSE(input_unit,STATUS='KEEP') ELSE WRITE(0,*)"Unable to OPEN file";STOP END IF END IF WRITE(0,*)'READ ',TRIM(namec),' and ',TRIM(namep) END IF DO jb=1,jbMax !pi --- 180 ! y x ! ! y = X*pi/180 !radianos ! lati (jbMax+1-jb) = (latic(1,jb)+90.0_r8)*pi/180.0_r8 xlati (jbMax+1-jb) = (latic(1,jb))*pi/180.0_r8 DO ib=1,nCols lonrad(ib) = longc(ib,1) gu (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*uvesc(ib,jb) & + ((step*dt)/(DELTAIN))*uvesp(ib,jb) gv (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*vvesc(ib,jb) & + ((step*dt)/(DELTAIN))*vvesp(ib,jb) gt (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) gq (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) & + ((step*dt)/(DELTAIN))*q02mp(ib,jb) gps (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*pslcc(ib,jb) & + ((step*dt)/(DELTAIN))*pslcp(ib,jb) var (ib, jb) = ((DELTAIN-step*dt)/(DELTAIN))*precc(ib,jb) & + ((step*dt)/(DELTAIN))*precp(ib,jb) dcupr_sib (ib,jbMax+1-jb) = (var(ib,jb)/86400.0_r8)*1.0e+3_r8 dlspr_sib (ib,jbMax+1-jb) = 0.0_r8 swdown (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*ocisc(ib,jb) & + ((step*dt)/(DELTAIN))*ocisp(ib,jb) dlwbot_sib(ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) & + ((step*dt)/(DELTAIN))*olisp(ib,jb) sigki = 1.0_r8 delsig = 0.010000_r8 tseam (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) tsea (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) slrad (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*olisc(ib,jb) & + ((step*dt)/(DELTAIN))*olisp(ib,jb) ! zorl........zorl (i)= 100.0 *zgrav*speedm(i)*rhi(i) ! zgrav =0.032 /grav zorl (ib,jbMax+1-jb) = 100.0_r8*(0.32_r8/9.8_r8)*(sqrt( gu(ib,jbMax+1-jb)**2 + gv(ib,jbMax+1-jb)**2))*0.01 tsurf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*t02mc(ib,jb) & + ((step*dt)/(DELTAIN))*t02mp(ib,jb) qsurf (ib,jbMax+1-jb) = ((DELTAIN-step*dt)/(DELTAIN))*q02mc(ib,jb) & + ((step*dt)/(DELTAIN))*q02mp(ib,jb) tmtx (ib,:,jbMax+1-jb) = 0.010000_r8 qmtx (ib,:,jbMax+1-jb) = 0.010000_r8 umtx (ib,:,jbMax+1-jb) = 0.010000_r8 END DO END DO step=step+1 ! ! mon is the month used for vegetation data input ! DO jb=1,jbMax !DO ib=1,ibMax,1 ! month(ib,jb)=idatec(2) ! IF((((lati(jb)*180.0_r8)/3.1415926e0_r8)-90.0_r8) > 0.0_r8 ) THEN ! month(ib,jb) = month(ib,jb) + 6 ! IF(month(ib,jb).GE.13) month(ib,jb) = month(ib,jb)-12 ! END IF !END DO ! ! computation of astronomical parameters ! sdelt ;solar inclination ! etime ;correction factor to local time ! ratio ;factor relating to the distance between the earth and the sun ! CALL radtim(idatec,sdelt ,ratio ,etime ,tod ,xday ,yrl) sw=1370.0_r8*ratio sindel = SIN(sdelt) cosdel = COS(sdelt) fimxi = 24.0e0_r8 /360.0_r8 ctime = alon/15.0e0_r8 cos2 = 0.0e0_r8 frh = ( MOD(tod+0.03125_r8,3600.0_r8)-0.03125_r8)/3600.0_r8 DO ib=1,ibMax zenith1(ib) = sindel*COS(lati(jb)) ENDDO DO ib=1,ibMax btime = fimxi*lonrad(ib)+ctime atime = etime+pai12*(12.0_r8-idatec(1)-frh-btime) zenith2 (ib) = cosdel*SIN(lati(jb))*COS(atime) coszen (ib,jb) = MAX (0.0_r8,zenith1(ib) + zenith2(ib)) sunang = MAX( 0.01_r8, coszen (ib,jb) ) cloud(ib,jb) = (sw * sunang - swdown(ib,jb)) / (sw * sunang) cloud(ib,jb) = MAX(cloud(ib,jb),0.0_r8) cloud(ib,jb) = MIN(cloud(ib,jb),1.0_r8) !cloud(ib,jb) = MAX(0.58_r8 ,cloud(ib,jb)) END DO END DO ! ! do for each waveband ! DO ib = 1, nband DO jb=1,jbMax DO i=1,nCols ! calculate the solar transmission through the atmosphere ! using simple linear function of tranmission and cloud cover ! ! note that the 'cloud cover' data is typically obtained from ! sunshine hours -- not direct cloud observations ! ! where, cloud cover = 1 - sunshine fraction ! ! different authors present different values for the slope and ! intercept terms of this equation ! ! Friend, A: Parameterization of a global daily weather generator for ! terrestrial ecosystem and biogeochemical modelling, Ecological ! Modelling ! ! Spitters et al., 1986: Separating the diffuse and direct component ! of global radiation and its implications for modeling canopy ! photosynthesis, Part I: Components of incoming radiation, ! Agricultural and Forest Meteorology, 38, 217-229. ! ! A. Friend : trans = 0.251 + 0.509 * (1.0 - cloud(i)) ! Spitters et al. : trans = 0.200 + 0.560 * (1.0 - cloud(i)) ! ! we are using the values from A. Friend ! trans = 0.251_r8 + 0.509_r8 * (1.0_r8 - cloud(i,jb)) ! ! calculate the fraction of indirect (diffuse) solar radiation ! based upon the cloud cover ! ! note that these relationships typically are measured for either ! monthly or daily timescales, and may not be exactly appropriate ! for hourly calculations -- however, in ibis, cloud cover is fixed ! through the entire day so it may not make much difference ! ! method i -- ! ! we use a simple empirical relationships from Nikolov and Zeller (1992) ! ! Nikolov, N. and K.F. Zeller, 1992: A solar radiation algorithm for ecosystem ! dynamics models, Ecological Modelling, 61, 149-168. ! fdiffuse = 1.0045_r8 + 0.0435_r8 * trans - 3.5227_r8 * trans**2 + 2.6313_r8 * trans**3 ! IF (trans.GT.0.75_r8) fdiffuse = 0.166_r8 ! ! method ii -- ! ! another method was suggested by Spitters et al. (1986) based on ! long-term data from the Netherlands ! ! Spitters et al., 1986: Separating the diffuse and direct component ! of global radiation and its implications for modeling canopy ! photosynthesis, Part I: Components of incoming radiation, ! Agricultural and Forest Meteorology, 38, 217-229. ! ! if ((trans.eq.0.00).and.(trans.lt.0.07)) then ! fdiffuse = 1.0 ! else if ((trans.ge.0.07).and.(trans.lt.0.35)) then ! fdiffuse = 1.0 - 2.3 * (trans - 0.07)**2 ! else if ((trans.ge.0.35).and.(trans.lt.0.75)) then ! fdiffuse = 1.33 - 1.46 * trans ! else ! fdiffuse = 0.23 ! endif ! ! do for each waveband ! ! ! calculate the fraction in each waveband ! ! frac = 0.46 + 0.08 * float(ib - 1) frac = 0.827_r8 + 0.146_r8 * REAL((ib - 1),KIND=r8) ! ! calculate the direct and indirect solar radiation ! solad(i,ib) = sw * coszen(i,jb) * frac * trans * (1.0_r8 - fdiffuse) ! solai(i,ib) = sw * coszen(i,jb) * frac * trans * fdiffuse cosz_sib (i,jb) =coszen (i,jb) radvbc_sib (i,jb) =solad (i,1) radnbc_sib (i,jb) =solad (i,2) radvdc_sib (i,jb) =solai (i,1) radndc_sib (i,jb) =solai (i,2) ! END DO END DO END DO RETURN 100 STOP END SUBROUTINE RunAtmosModel ! radtim :calculates the astronomical parameters: solar inclination, ! correction factor to local time, factor relating to the distance ! between earth and sun, and the julian day. SUBROUTINE radtim (id ,delta ,ratio ,etime ,tod ,xday ,yrl) ! !========================================================================== ! !========================================================================== ! id(1)....hour(00/12) ! id(2)....month ! id(3)....day of month ! id(4)....year ! delta....solar inclination ! ratio....factor relating to the distance between the earth and the sun ! etime....correction factor to local time ! tod......model forecast time of day in seconds ! xday.....is julian day - 1 with fraction of day ! pai......constant pi=3.1415926 ! yrl......length of year in days ! monl.....length of each month in days !========================================================================== INTEGER, INTENT(in ) :: id(4) REAL(KIND=r8), INTENT(out) :: delta REAL(KIND=r8), INTENT(out) :: ratio REAL(KIND=r8), INTENT(out) :: etime REAL(KIND=r8), INTENT(in ) :: tod REAL(KIND=r8), INTENT(out) :: xday REAL(KIND=r8), INTENT(in ) :: yrl REAL(KIND=r8), PARAMETER :: day0=-1.0_r8 REAL(KIND=r8), PARAMETER :: f3600=3.6e3_r8 REAL(KIND=r8) :: psi ! ! id is now assumed to be the current date and hour ! xday=id(1)*f3600 xday=xday+MOD(tod,f3600) xday=monday(id(2))+id(3)+xday/86400.0_r8 IF (yrl == 365.25e0_r8) THEN xday=xday-MOD(id(4)+3,4)*0.25_r8 IF(MOD(id(4),4).EQ.0.AND.id(2).GT.2)xday=xday+1.0e0_r8 END IF xday= MOD(xday-1.0e0_r8,yrl) IF (xday > day0) THEN psi=2.0e0_r8*pai*xday/yrl delta=0.006918e0_r8-0.399912e0_r8*COS( psi)+0.070257e0_r8*SIN(psi) & -0.006758e0_r8*COS(2.0e0_r8*psi)+0.000907e0_r8*SIN(2.0e0_r8*psi) & -0.002697e0_r8*COS(3.0e0_r8*psi)+0.001480e0_r8*SIN(3.0e0_r8*psi) ratio=1.000110e0_r8+0.034221e0_r8*COS( psi)+0.001280e0_r8*SIN(psi) & +0.000719e0_r8*COS(2.0e0_r8*psi)+0.000077e0_r8*SIN(2.0e0_r8*psi) etime=0.000075e0_r8+0.001868e0_r8*COS( psi)-0.032077e0_r8*SIN(psi) & -0.014615e0_r8*COS(2.0e0_r8*psi)-0.040849e0_r8*SIN(2.0e0_r8*psi) ELSE WRITE(0,20)id,tod,xday WRITE(0,20)id,tod,xday STOP 2020 END IF 20 FORMAT(' BAD DATE IN RADTIM. ID=',4I5,' TOD=',G16.8,' XDAY=', & G16.8) END SUBROUTINE radtim SUBROUTINE InitRadtim(monl) INTEGER, INTENT(in ) :: monl(12) INTEGER :: m monday(1)=0 DO m=2,12 monday(m)=monday(m-1)+monl(m-1) END DO END SUBROUTINE InitRadtim END MODULE AtmosModel