!WRF:MODEL_RA:RADIATION
!
MODULE module_ra_gfdleta 2
INTEGER, PARAMETER :: NL=81
INTEGER, PARAMETER :: NBLY=15
REAL , SAVE, DIMENSION(37,NL) :: XDUO3N,XDO3N2,XDO3N3,XDO3N4
REAL , SAVE, DIMENSION(NL) :: PRGFDL
REAL , SAVE :: AB15WD,SKO2D,SKC1R,SKO3R
REAL , SAVE :: EM1(28,180),EM1WDE(28,180),TABLE1(28,180), &
TABLE2(28,180),TABLE3(28,180),EM3(28,180), &
SOURCE(28,NBLY), DSRCE(28,NBLY)
REAL ,SAVE, DIMENSION(5040):: T1,T2,T4,EM1V,EM1VW,EM3V
REAL ,SAVE :: R1
! Created by CO2 initialization
REAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: CO251,CDT51,CDT58,C2D51,&
C2D58,CO258
REAL, SAVE, ALLOCATABLE, DIMENSION(:) :: STEMP,GTEMP,CO231,CO238, &
C2D31,C2D38,CDT31,CDT38, &
CO271,CO278,C2D71,C2D78, &
CDT71,CDT78
REAL, SAVE, ALLOCATABLE, DIMENSION(:) :: CO2M51,CO2M58,CDTM51,CDTM58, &
C2DM51,C2DM58
! Used by CO2 initialization
! COMMON/PRESS/PA(109)
! COMMON/TRAN/ TRANSA(109,109)
! COMMON/COEFS/XA(109),CA(109),ETA(109),SEXPV(109),CORE,UEXP,SEXP
REAL ,SAVE, DIMENSION(109) :: PA, XA, CA, ETA, SEXPV
REAL ,SAVE, DIMENSION(109,109) :: TRANSA
REAL ,SAVE :: CORE,UEXP,SEXP
#ifndef TRIEDNTRUE
! EQUIVALENCE (EM1V(1),EM1(1,1)),(EM1VW(1),EM1WDE(1,1))
! EQUIVALENCE (EM3V(1),EM3(1,1))
! EQUIVALENCE (T1(1),TABLE1(1,1)),(T2(1),TABLE2(1,1)), &
! (T4(1),TABLE3(1,1))
#endif
CONTAINS
!---------------------------------------------------------------------
SUBROUTINE GFDLETAINIT(SFULL,SHALF,PPTOP,JULYR,JULDAY,GMT, & 2,4
kds, kde, kms, kme, kts, kte)
IMPLICIT NONE
INTEGER, INTENT(IN) :: kds, kde, kms, kme, kts, kte
REAL, DIMENSION(kms:kme), INTENT(IN) :: SFULL, SHALF
INTEGER, INTENT(IN) :: JULDAY,JULYR
REAL, INTENT(IN) :: GMT, PPTOP
#ifndef TRIEDNTRUE
INTEGER IHRST
CALL CO2O3
(sfull,shalf,pptop,kme-kms,kme-kms+1,kme-kms+2)
CALL O3CLIM
CALL TABLE
IHRST = NINT(GMT)
CALL SOLARD(IHRST,JULDAY,JULYR)
#endif
END SUBROUTINE GFDLETAINIT
!---------------------------------------------------------------------
SUBROUTINE ETARA (DT,THRATEN,THRATENLW,THRATENSW,PI3D & 2,2
,XLAND,p8w,dz8w,RHO_PHY,P_PHY,T &
,QV,QL,TSK2D,GLW,GSW &
,GLAT,GLON,HTOP,HBOT,ALBEDO,CUPPT &
,VEGFRA,SNOW,G,GMT,NSTEPRA,itimestep &
,julyr,julday,gfdl_lw,gfdl_sw &
,IDS,IDE,JDS,JDE,KDS,KDE &
,IMS,IME,JMS,JME,KMS,KME &
,ITS,ITE,JTS,JTE,KTS,KTE )
!---------------------------------------------------------------------
IMPLICIT NONE
!---------------------------------------------------------------------
INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
,IMS,IME,JMS,JME,KMS,KME &
,ITS,ITE,JTS,JTE,KTS,KTE,itimestep,NSTEPRA
INTEGER,INTENT(IN) :: julyr,julday
REAL ,INTENT(IN) :: DT,GMT,G
REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme):: &
THRATEN,THRATENLW,THRATENSW
REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme)::p8w,dz8w, &
rho_phy,p_phy,PI3D
REAL, INTENT(IN), DIMENSION(ims:ime, jms:jme):: XLAND, TSK2D, &
ALBEDO,VEGFRA,SNOW
REAL, INTENT(IN), DIMENSION(its:ite, jts:jte):: GLAT,GLON
REAL, INTENT(INOUT), DIMENSION(ims:ime, jms:jme):: HTOP,HBOT,CUPPT
REAL, DIMENSION(ims:ime, jms:jme):: GLW,GSW
REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme):: QV,QL,T
! REAL, INTENT(IN), DIMENSION(37*kte) :: RAD1,RAD2,RAD3,RAD4
LOGICAL, INTENT(IN) :: gfdl_lw,gfdl_sw
REAL, DIMENSION(its:ite, kms:kme, jts:jte)::TFLIP,QFLIP,QLFLIP,PFLIP
REAL, DIMENSION(its:ite, kms:kme, jts:jte)::P8WFLIP,PHYD
REAL, DIMENSION(its:ite, kts:kte, jts:jte)::TENDS,TENDL
INTEGER :: IDAT(3),Jmonth,Jday
INTEGER :: I,J,K,KFLIP,IHRST
#ifndef TRIEDNTRUE
IF (gfdl_lw .and. gfdl_sw ) goto 100
! NEED HYDROSTATIC PRESSURE HERE (MONOTONIC CHANGE WITH HEIGHT)
DO J = jts,jte
DO I = its,ite
PHYD(I,KTS,J)=P8W(I,KTS,J)
ENDDO
ENDDO
DO K = KTS,KTE
DO J = jts,jte
DO I = its,ite
PHYD(I,K+1,J) = PHYD(I,K,J) - G*RHO_PHY(I,K,J)*DZ8W(I,K,J)
ENDDO
ENDDO
ENDDO
DO K = KMS,KME
KFLIP=KME+1-K
DO J = jts,jte
DO I = its,ite
P8WFLIP(I,K,J) = PHYD(I,KFLIP,J)
ENDDO
ENDDO
ENDDO
DO K = KTS,KTE
KFLIP=KTE+1-K
DO J = jts,jte
DO I = its,ite
TFLIP (I,K,J) = T(I,KFLIP,J)
QFLIP (I,K,J) = MAX(0.,QV(I,KFLIP,J)/(1.+QV(I,KFLIP,J)))
QLFLIP(I,K,J) = QL(I,KFLIP,J)
! PFLIP (I,K,J) = P_PHY(I,KFLIP,J)
! USE MONOTONIC HYDROSTATIC PRESSURE INTERPOLATED TO MID-LEVEL
PFLIP (I,K,J) = 0.5*(P8WFLIP(I,K,J)+P8WFLIP(I,K+1,J))
ENDDO
ENDDO
ENDDO
CALL CAL_MON_DAY(JULDAY,julyr,Jmonth,Jday)
IDAT(1) = Jmonth
IDAT(2) = Jday
IDAT(3) = julyr
IHRST = nint(GMT)
! CALL SOLARD(R1,IHRST,IDAT)
! CALL SOLARD(R1,IHRST,julday)
CALL RADTN (DT,TFLIP,QFLIP,QLFLIP,PFLIP,P8WFLIP,XLAND,TSK2D, &
GLAT,GLON,HTOP,HBOT,ALBEDO,CUPPT, &
! IHE,IHW,ACFRCV,NCFRCV,ACFRST,NCFRST,VEGFRC, &
VEGFRA,SNOW,GLW,GSW, &
! IDAT,IHRST, &
IDAT,IHRST, &
! R1,NRADS,NRADL,RESTRT,NTSD,NPHS, &
NSTEPRA,NSTEPRA,itimestep,1, &
TENDS,TENDL, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
IF ( gfdl_lw ) then
DO J=JTS,JTE
DO K = KTS,KTE
KFLIP=KTE+1-K
DO I=ITS,ITE
THRATENLW(I,K,J)=TENDL(I,KFLIP,J)/PI3D(I,K,J)
THRATENSW(I,K,J)=TENDS(I,KFLIP,J)/PI3D(I,K,J)
THRATEN(I,K,J) =THRATEN(I,K,J) + THRATENLW(I,K,J)
ENDDO
ENDDO
ENDDO
ENDIF
! This assumed that longwave is called first in the radition_driver.
IF ( gfdl_sw ) then
DO J=JTS,JTE
DO K = KTS,KTE
KFLIP=KTE+1-K
DO I=ITS,ITE
THRATENSW(I,K,J)=TENDS(I,KFLIP,J)/PI3D(I,K,J)
ENDDO
ENDDO
ENDDO
ENDIF
100 IF ( gfdl_sw ) then
DO J=JTS,JTE
DO K = KTS,KTE
KFLIP=KTE+1-K
DO I=ITS,ITE
THRATEN(I,K,J) =THRATEN(I,K,J) + THRATENSW(I,K,J)
ENDDO
ENDDO
ENDDO
ENDIF
#endif
END SUBROUTINE ETARA
#ifndef TRIEDNTRUE
!----------------------------------------------------------------------
SUBROUTINE RADTN(DT,T,Q,CWM,PFLIP,P8WFLIP,XLAND,TSK2D, & 1,10
GLAT,GLON,HTOP,HBOT,ALB,CUPPT, &
! RAD1,RAD2,RAD3,RAD4, &
! TABLE1,TABLE2,TABLE3,EM1,EM1WDE,EM3, &
VEGFRC,SNO,GLW,GSW, &
! IDAT,LTOP,IHRST,PRGFDL, &
IDAT,IHRST, &
NRADS,NRADL,NTSD,NPHS, &
! SKO3R,AB15WD,SKC1R,SKO2D, &
TENDS,TENDL, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!----------------------------------------------------------------------
IMPLICIT NONE
!----------------------------------------------------------------------
! GLAT : geodetic latitude in radians of the mass points on the computational grid.
! CZEN : instantaneous cosine of the solar zenith angle.
! HTOP : (REAL) model layer number that is highest in the atmosphere
! in which convective cloud occurred since the previous call to the
! radiation driver.
! HBOT : (REAL) model layer number that is lowest in the atmosphere
! in which convective cloud occurred since the previous call to the
! radiation driver.
! ALB : is no longer used in the operational radiation. Prior to 24 July,
! ALB was the climatological albedo that was modified within RADTN to
! account for vegetation fraction and snow.
! CUPPT: accumulated convective precipitation (meters) since the
! last call to the radiation.
! THS : potential temperature of the ground surface.
! IHE and IHW are relative location indices needed to locate neighboring
! points on the Eta's Arakawa E grid since arrays are indexed locally on
! each MPI task rather than globally. IHE refers to the adjacent grid
! point (a V point) to the east of the mass point being considered. IHW
! is the adjacent grid point to the west of the given mass point.
! IRAD is a relic from older code that is no longer needed.
! ACFRCV : sum of the convective cloud fractions that were computed
! during each call to the radiation between calls to the subroutines that
! do the forecast output.
! NCFRCV : the total number of times in which the convective cloud
! fraction was computed to be greater than zero in the radiation between
! calls to the output routines. In the post-processor, ACFRCV is divided
! by NCFRCV to yield an average convective cloud fraction.
! ACFRST and NCFRST are the analogs for stratiform cloud cover.
! VEGFRC is the fraction of the gridbox with vegetation.
! LVL holds the number of model layers that lie below the ground surface
! at each point. Clearly for sigma coordinates LVL is zero everywhere.
! CTHK : an assumed maximum thickness of stratiform clouds currently set
! to 20000 Pascals. I think this is relevant for computing "low",
! "middle", and "high" cloud fractions which are post-processed but which
! do not feed back into the integration.
! IDAT : a 3-element integer array holding the month, day, and year,
! respectively, of the date for the start time of the free forecast.
! ABCFF : holds coefficients for various absorption bands. You can see
! where they are set in GFDLRD.F.
! LTOP : a 3-element integer array holding the model layer that is at or
! immediately below the specified pressure levels for the tops
! of "high" (15000 Pa), "middle" (35000 Pa), and "low" (64200 Pa)
! stratiform clouds. These are for the diagnostic cloud layers
! needed in the output but not in the integration.
! R1 : earth-sun distance in astronomical units.
! NRADS : integer number of fundamental timesteps (our smallest
! timestep, i.e., the one for inertial gravity wave adjustment)
! between updates of the shortwave tendencies. Currently we
! update the shortwave every hour.
! NRADL : integer number of fundamental timesteps between updates of
! the longwave tendencies. Currently we update the longwave
! every two hours.
! NTSD : integer counter of the fundamental timesteps that have
! elapsed since the start of the forecast.
!-----------------------------------------------------------------
! INTEGER, PARAMETER :: NL=81
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
INTEGER, INTENT(IN) :: NRADS,NRADL,NTSD,NPHS
! LOGICAL, INTENT(IN) :: RESTRT
REAL , INTENT(IN) :: DT
! REAL , INTENT(IN), DIMENSION(37,NL) :: XDUO3N,XDO3N2,XDO3N3,XDO3N4
INTEGER, INTENT(IN), DIMENSION(3) :: IDAT
!----------------------------------------------------------------------
REAL, PARAMETER :: CAPA=0.28589641,RTD=57.2957795, WA=.10, &
WG=1.-WA,KSMUD=0
REAL, PARAMETER :: A1=610.78,A2=17.2693882,A3=273.16, &
A4=35.86, PQ0=379.90516,SNOALB=0.55
INTEGER :: LM1,LP1,LM
INTEGER, INTENT(IN) :: IHRST
! REAL, INTENT(IN), DIMENSION(NL) :: PRGFDL
!
REAL, PARAMETER :: SLPM=1.01325E5,EPSQ1=1.E-5,EPSQ=2.E-12, &
EPSO3=1.E-10,HPINC=1.E1, CLDRH0=0.80, &
TRESH=1.00,RNRM=1./(TRESH-CLDRH0), &
CLDRH2=0.90,TRESH2=1.00, &
RNRM2=1./(TRESH2-CLDRH2), CLAPSE=-0.0005, &
CLPSE=-0.0006,DCLPS=-0.0001, CM1=2937.4, &
CM2=4.9283,CM3=23.5518,EPS=0.622, &
PBOT=10000.0, STBOL=5.67E-8, &
PI2=2.*3.14159265,RLAG=14.8125
!
INTEGER, PARAMETER :: NB=12
INTEGER, PARAMETER :: K15=SELECTED_REAL_KIND(15)
REAL (KIND=K15) :: DDX,EEX,PROD
! REAL :: PROD,DDX,EEX
! REAL, INTENT(IN) :: SKO3R,AB15WD,SKC1R,SKO2D
!-----------------------------------------------------------------------
LOGICAL :: SHORT,LONG
LOGICAL :: BCLD(its:ite),BTEMP1(its:ite)
LOGICAL :: BITX,BITY,BITZ,BITW,BIT1,BIT2,BITC,BITS,BITCP1,BITSP1
LOGICAL :: CNCLD
!-----------------------------------------------------------------------
REAL, INTENT(IN), DIMENSION(ims:ime,jms:jme) :: XLAND,TSK2D
REAL, INTENT(IN), DIMENSION(its:ite, kms:kme, jts:jte):: Q,CWM,T
REAL, INTENT(IN), DIMENSION(its:ite, kms:kme, jts:jte):: PFLIP,P8WFLIP
! REAL, INTENT(IN), DIMENSION(28,180) :: TABLE1,TABLE2,TABLE3,EM3,EM1,EM1WDE
REAL, INTENT(OUT), DIMENSION(ims:ime, jms:jme):: GLW,GSW
! REAL, INTENT(IN), DIMENSION(kms:kme) :: ETAD
! REAL, INTENT(IN), DIMENSION(kms:kme) :: AETA
!-----------------------------------------------------------------------
REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: HTOP,HBOT
REAL, INTENT(IN ), DIMENSION(ims:ime,jms:jme) :: ALB,SNO
REAL, INTENT(IN ), DIMENSION(its:ite,jts:jte) :: GLAT,GLON
!-----------------------------------------------------------------------
REAL, DIMENSION(its:ite,jts:jte) :: CZEN
REAL, DIMENSION(its:ite,jts:jte) :: CZMEAN,SIGT4
INTEGER, DIMENSION(its:ite, jts:jte):: LMH
REAL, DIMENSION(its:ite, jts:jte):: U00
!-----------------------------------------------------------------------
REAL, DIMENSION(2*kte) :: UL
!-----------------------------------------------------------------------
REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: CUPPT
REAL, DIMENSION(its:ite,jts:jte) :: CFRACL,CFRACM,CFRACH
! REAL, DIMENSION(37*kte) :: RAD1,RAD2,RAD3,RAD4
!-----------------------------------------------------------------------
! INTEGER,INTENT(IN), DIMENSION(jms:jme) :: IHE,IHW
!-----------------------------------------------------------------------
! REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: ACFRCV,ACFRST
! INTEGER, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: NCFRCV,NCFRST
!----------------------------------------------------------------------
REAL, DIMENSION(its:ite,jts:jte) :: RLWIN,RLWOUT,RLWTOA
!----------------------------------------------------------------------
REAL, INTENT(IN), DIMENSION(ims:ime,jms:jme) :: VEGFRC
REAL, INTENT(INOUT),DIMENSION(its:ite,kts:kte,jts:jte) :: TENDS,TENDL
REAL, DIMENSION(its:ite,jts:jte) :: RSWIN,RSWOUT,RSWTOA
REAL, DIMENSION(its:ite,kts:kte,jts:jte):: RSWTT,RLWTT
!-----------------------------------------------------------------------
REAL :: CTHK(3)
DATA CTHK/20000.0,20000.0,20000.0/
! INTEGER,INTENT(IN) :: LTOP(3)
REAL :: PTOPC(4)
REAL, DIMENSION(9) :: CC,PPT
!-----------------------------------------------------------------------
REAL :: ABCFF(NB)
INTEGER,DIMENSION(its:ite,jts:jte) :: LVL
REAL, DIMENSION(its:ite, jts:jte):: PDSL,FNE,FSE,TL
REAL, DIMENSION( 0:kte) :: CLDAMT
REAL, DIMENSION(its:ite,3):: CLDCFR
INTEGER, DIMENSION(its:ite,3):: MBOT,MTOP
REAL, DIMENSION(its:ite) :: PSFC,TSKN,ALBEDO,XLAT,COSZ, &
SLMSK,CV,SV,FLWUP, &
FSWDN,FSWUP,FSWDNS,FSWUPS,FLWDNS, &
FLWUPS
INTEGER, DIMENSION(its:ite) :: ICVB,ICVT
REAL, DIMENSION(its:ite,kts:kte) :: PMID,TMID
REAL, DIMENSION(its:ite,kts:kte) :: QMID,THMID,OZN,POZN
REAL, DIMENSION(its:ite,jts:jte) :: TOT
REAL, DIMENSION(its:ite,kts:kte+1) :: PINT,EMIS,CAMT
INTEGER,DIMENSION(its:ite,kts:kte+1) :: ITYP,KBTM,KTOP
INTEGER,DIMENSION(its:ite) :: NCLDS,KCLD
REAL, DIMENSION(its:ite) :: CSTR,TAUC,TAUDAR
REAL, DIMENSION(its:ite,NB,kts:kte+1) ::RRCL,TTCL
INTEGER, DIMENSION(its:ite,kts:kte):: IW
REAL, DIMENSION(its:ite,kts:kte):: CCR,CSMID,WMID,HMID,CCMID
REAL, DIMENSION(its:ite) :: BMID,UMID
REAL :: PLOMD,PMDHI,PHITP,P400,PLBTM
INTEGER :: NFILE
REAL :: UTIM,CLSTP,TIME,DAYI,HOUR,ADDL,RANG,RSIN1,RCOS1,RCOS2
REAL :: TIMES,EXNER,APES,SNOFAC,US,CCLIMIT,CLIMIT,HH,TKL,QKL
REAL :: CWMKL,TMT15,AI,BI,PP,QW,P1,CC2,CC1,PMOD,CLPFIL,DTHDP,DDP
REAL :: ARG,RQKL,FIW,U00KL,QI,FIQ,TMT0,P2,QINT,QC,CLDMAX
REAL :: CL1,CL2,CR1,DPCL,QSUM,PRS1,PRS2,DELP,TCLD,DD,EE,AA,FF
REAL :: BB,GG,DENOM,FCTRA,FCTRB,PDSLIJ,CFRAVG,SNOMM
INTEGER :: I,J,MYJS,MYJE,MYIS,MYIE,NTSPH,NRADPP,ITIMSW,ITIMLW,JD,II
INTEGER :: L,N,LML,LVLIJ,IR,KNTLYR,LL,NC,L400,NMOD,LTROP,IWKL
INTEGER :: NLVL,MALVL,LLTOP,LLBOT,KBT2,KTH1,KBT1,KTH2,KTOP1
INTEGER :: NBAND,NCLD,LBASE,NKTP,NBTM,KS,MYJS1,MYJS2,MYJE2,MYJE1
! REAL,DIMENSION(5040):: T1,T2,T4,EM1V,EM1VW,EM3V
! EQUIVALENCE (EM1V(1),EM1(1,1)),(EM1VW(1),EM1WDE(1,1))
! EQUIVALENCE (EM3V(1),EM3(1,1))
! EQUIVALENCE (T1(1),TABLE1(1,1)),(T2(1),TABLE2(1,1)), &
! (T4(1),TABLE3(1,1))
DATA PLOMD/64200./,PMDHI/35000./,PHITP/15000./,P400/40000./, &
PLBTM/105000./
DATA NFILE/14/
DATA CC/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8/
DATA PPT/.14,.31,.70,1.6,3.4,7.7,17.,38.,85./
DATA ABCFF/2*4.0E-5,.002,.035,.377,1.95,9.40,44.6,190.,989.,2706.,39011./
!----------------------------------------------------------------------
MYJS=jts
MYJE=jte
MYIS=its
MYIE=ite
MYJS1=jts !????
MYJE1=jte
MYJS2=jts
MYJE2=jte
LM=kte
LM1=LM-1
LP1=LM+1
UTIM=1.
CNCLD=.TRUE.
do j = jts,jte
do i = its,ite
LMH(i,j) = kme-1
LVL(i,j) = 0
enddo
enddo
do j = jts,jte
do i = its,ite
U00(I,J)=(2.-XLAND(i,j))*0.75+(XLAND(i,j) -1.)*0.80
enddo
enddo
DO L=1,2*LM
IF(L.GE.LM-10.AND.L.LE.LM)THEN
UL(L)=0.1*FLOAT(L-LM+10)
ELSE
UL(L)=0.
ENDIF
ENDDO
!***
!*** ASSIGN THE PRESSURES FOR CLOUD DOMAIN BOUNDARIES
!***
PTOPC(1)=PLBTM
PTOPC(2)=PLOMD
PTOPC(3)=PMDHI
PTOPC(4)=PHITP
!***
!*** FIND THE 'SEA LEVEL PRESSURE'.
!***
! DO J=MYJS,MYJE
! DO I=MYIS,MYIE
! PDSL(I,J)=PD(I,J)
! ENDDO
! ENDDO
!**********************************************************************
!*** THE FOLLOWING CODE IS EXECUTED EACH TIME THE RADIATION IS CALLED.
!**********************************************************************
!----------------------CONVECTION--------------------------------------
! NRADPP IS THE NUMBER OF TIME STEPS TO ACCUMULATE CONVECTIVE PRECIP
! FOR RADIATION
! NOTE: THIS WILL NOT WORK IF NRADS AND NRADL ARE DIFFERENT UNLESS
! THEY ARE INTEGER MULTIPLES OF EACH OTHER
! CLSTP IS THE NUMBER OF HOURS OF THE ACCUMULATION PERIOD
!
NTSPH=NINT(3600./DT)
NRADPP=MIN(NRADS,NRADL)
CLSTP=1.0*NRADPP/NTSPH
!----------------------CONVECTION--------------------------------------
!***
!*** STATE WHETHER THE SHORT OR LONGWAVE COMPUTATIONS ARE TO BE DONE.
!***
!SH SHORT=.FALSE.
!SH LONG=.FALSE.
!SH IF(MOD(NTSD,NRADS).EQ.1.OR.RESTRT)SHORT=.TRUE.
!SH IF(MOD(NTSD,NRADL).EQ.1.OR.RESTRT)LONG=.TRUE.
!SH IF(NTSD .eq. 1 .or. MOD((NTSD),NRADS).EQ.0.)SHORT=.TRUE.
!SH IF(NTSD .eq. 1 .or. MOD((NTSD),NRADL).EQ.0.)LONG=.TRUE.
SHORT=.TRUE.
LONG=.TRUE.
ITIMSW=0
ITIMLW=0
IF(SHORT)ITIMSW=1
IF(LONG) ITIMLW=1
!***
!*** FIND THE MEAN COSINE OF THE SOLAR ZENITH ANGLE
!*** BETWEEN THE CURRENT TIME AND THE NEXT TIME RADIATION IS
!*** CALLED. ONLY AVERAGE IF THE SUN IS ABOVE THE HORIZON.
!***
TIME=(NTSD-1)*DT
! CALL ZENITH(TIME,DAYI,HOUR)
CALL ZENITH(TIME,DAYI,HOUR,IDAT,IHRST,GLON,GLAT,CZEN, &
MYIS,MYIE,MYJS,MYJE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
JD=INT(DAYI+0.50)
ADDL=0.
IF(MOD(IDAT(3),4).EQ.0)ADDL=1.
RANG=PI2*(DAYI-RLAG)/(365.25+ADDL)
RSIN1=SIN(RANG)
RCOS1=COS(RANG)
RCOS2=COS(2.*RANG)
IF(SHORT)THEN
DO J=MYJS,MYJE
DO I=MYIS,MYIE
CZMEAN(I,J)=0.
TOT(I,J)=0.
ENDDO
ENDDO
!
DO II=0,NRADS,NPHS
TIMES=(NTSD-1)*DT+II*DT
! CALL ZENITH(TIMES,DAYI,HOUR)
CALL ZENITH(TIMES,DAYI,HOUR,IDAT,IHRST,GLON,GLAT,CZEN, &
MYIS,MYIE,MYJS,MYJE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
DO J=MYJS,MYJE
DO I=MYIS,MYIE
IF(CZEN(I,J).GT.0.)THEN
CZMEAN(I,J)=CZMEAN(I,J)+CZEN(I,J)
TOT(I,J)=TOT(I,J)+1.
ENDIF
ENDDO
ENDDO
ENDDO
DO J=MYJS,MYJE
DO I=MYIS,MYIE
IF(TOT(I,J).GT.0.)CZMEAN(I,J)=CZMEAN(I,J)/TOT(I,J)
ENDDO
ENDDO
ENDIF
!
!********************************************************************
!*** THIS IS THE BEGINNING OF THE PRIMARY LOOP THROUGH THE DOMAIN
!********************************************************************
! *********************
DO 700 J = MYJS, MYJE
! *********************
!
DO 125 L=1,LM
DO I=MYIS,MYIE
! IR=IRAD(I)
TMID(I,L)=T(I,1,J)
QMID(I,L)=EPSQ
CSMID(I,L)=0.
WMID(I,L)=0.
CCMID(I,L)=0.
IW(I,L)=0.
CCR(I,L)=0.
HMID(I,L)=0.
OZN(I,L)=EPSO3
TENDS(I,L,J)=0.
TENDL(I,L,J)=0.
ENDDO
125 CONTINUE
!
DO 140 N=1,3
DO I=MYIS,MYIE
CLDCFR(I,N)=0.
MTOP(I,N)=0
MBOT(I,N)=0
ENDDO
140 CONTINUE
!***
!*** FILL IN WORKING ARRAYS WHERE VALUES AT L=LM ARE THOSE THAT
!*** ARE ACTUALLY AT ETA LEVEL L=LMH.
!***
DO 200 I=MYIS,MYIE
! IR=IRAD(I)
LML=LMH(I,J)
LVLIJ=LVL(I,J)
!tlb BMID(I)=HBM2(IR,J)
! BMID(I)=HBM2(I,J)
UMID(I)=U00(I,J)
!
DO L=1,LML
! PMID(I,L+LVLIJ)=AETA(L)*PDSL(I,J)+PT
PMID(I,L+LVLIJ)=PFLIP(I,L,J)
! PINT(I,L+LVLIJ+1)=ETAD(L+1)*PDSL(I,J)+PT
PINT(I,L+LVLIJ+1)=P8WFLIP(I,L+1,J)
EXNER=(1.E5/PMID(I,L+LVLIJ))**CAPA
TMID(I,L+LVLIJ)=T(I,L,J)
THMID(I,L+LVLIJ)=T(I,L,J)*EXNER
QMID(I,L+LVLIJ)=Q(I,L,J)
WMID(I,L+LVLIJ)=CWM(I,L,J)
HMID(I,L+LVLIJ)=1.
ENDDO
!***
!*** FILL IN ARTIFICIAL VALUES ABOVE THE TOP OF THE DOMAIN.
!*** PRESSURE DEPTHS OF THESE LAYERS IS 1 HPA.
!*** TEMPERATURES ABOVE ARE ALREADY ISOTHERMAL WITH (TRUE) LAYER 1.
!***
IF(LVLIJ.GT.0)THEN
KNTLYR=0
!
DO L=LVLIJ,1,-1
KNTLYR=KNTLYR+1
! PMID(I,L)=PT-REAL(2*KNTLYR-1)*0.5*HPINC
PMID(I,L)=P8WFLIP(I,1,J)-REAL(2*KNTLYR-1)*0.5*HPINC
PINT(I,L+1)=PMID(I,L)+0.5*HPINC
EXNER=(1.E5/PMID(I,L))**CAPA
THMID(I,L)=TMID(I,L)*EXNER
ENDDO
ENDIF
!
IF(LVLIJ.EQ.0) THEN
! PINT(I,1)=PT
PINT(I,1)=P8WFLIP(I,1,J)
ELSE
PINT(I,1)=PMID(I,1)-0.5*HPINC
ENDIF
200 CONTINUE
!***
!*** FILL IN THE SURFACE PRESSURE, SKIN TEMPERATURE, GEODETIC LATITUDE,
!*** ZENITH ANGLE, SEA MASK, AND ALBEDO. THE SKIN TEMPERATURE IS
!*** NEGATIVE OVER WATER.
!***
DO 250 I=MYIS,MYIE
! PSFC(I)=PD(I,J)+PT
PSFC(I)=P8WFLIP(I,kme,J)
APES=(PSFC(I)*1.E-5)**CAPA
! TSKN(I)=THS(I,J)*APES*(1.-2.*SM(I,J))
IF ((XLAND(I,J)-1.5) .gt. 0.) then
TSKN(I)=-TSK2D(I,J)
ELSE
TSKN(I)=TSK2D(I,J)
ENDIF
! TSKN(I)=THS(I,J)*APES*(1.-2.*(XLAND(I,J)-1.))
! SLMSK(I)=SM(I,J)
SLMSK(I)=XLAND(I,J)-1.
!
! SNO(I,J)=AMAX1(SNO(I,J),0.)
SNOMM=AMAX1(SNO(I,J),0.)
SNOFAC=AMIN1(SNOMM/0.02, 1.0)
ALBEDO(I)=ALB(I,J)+(1.0-0.01*VEGFRC(I,J))*SNOFAC*(SNOALB-ALB(I,J))
!
XLAT(I)=GLAT(I,J)*RTD
COSZ(I)=CZMEAN(I,J)
250 CONTINUE
!-----------------------------------------------------------------------
!*******************STRATIFORM CLOUD SECTION***************************
!-----------------------------------------------------------------------
! CALCULATE STRATIFORM CLOUD COVERAGE AT EACH MODEL GRID POINT WHICH
! WILL BE USED IN THE MODEL RADIATION PARAMETERIZATION SCHEME.
!-----------------------------------------------------------------------
US=1.
CCLIMIT=1.0E-3
CLIMIT =1.0E-20
!------------------QW, QI AND QINT--------------------------------------
DO 280 I=MYIS,MYIE
LML=LMH(I,J)
LVLIJ=LVL(I,J)
!
DO 275 L=1,LML
LL=L+LVLIJ
! HH=HMID(I,LL)*BMID(I)
HH=HMID(I,LL)
TKL=TMID(I,LL)
QKL=QMID(I,LL)
CWMKL=WMID(I,LL)
TMT0=(TKL-273.16)*HH
TMT15=AMIN1(TMT0,-15.)*HH
AI=0.008855
BI=1.
!
IF(TMT0.LT.-20.)THEN
AI=0.007225
BI=0.9674
ENDIF
!
PP=PMID(I,LL)
QW=HH*PQ0/PP*EXP(HH*A2*(TKL-A3)/(TKL-A4))
QI=QW*(BI+AI*AMIN1(TMT0,0.))
QINT=QW*(1.-0.00032*TMT15*(TMT15+15.))
IF(TMT0.LE.-40.) QINT=QI
!-------------------ICE-WATER ID NUMBER IW------------------------------
U00KL=UMID(I)+UL(L)*(0.95-UMID(I))*UTIM
IF(TMT0.LT.-15.0)THEN
FIQ=QKL-U00KL*QI
IF(FIQ.GT.0..OR.CWMKL.GT.CLIMIT)THEN
IW(I,LL)=1
ELSE
IW(I,LL)=0
ENDIF
ENDIF
!
IF(TMT0.GE.0.)THEN
IW(I,LL)=0
ENDIF
!
IF(TMT0.LT.0.0.AND.TMT0.GE.-15.0)THEN
IW(I,LL)=0
IF(IW(I,LL-1).EQ.1.AND.CWMKL.GT.CLIMIT) IW(I,LL)=1
ENDIF
!
IWKL=IW(I,LL)
!
!----------------THE SATURATION SPECIFIC HUMIDITY------------------------
!
FIW=FLOAT(IWKL)
QC=(1.-FIW)*QINT+FIW*QI
!
!----------------THE RELATIVE HUMIDITY----------------------------------
!
IF(QC.LE.EPSQ1.OR.QKL.LE.EPSQ1)THEN
RQKL=0.
ELSE
RQKL=QKL/QC
ENDIF
!
!----------------CLOUD COVER RATIO CCR----------------------------------
!
IF(RQKL.GE.0.9999)THEN
CCR(I,LL)=AMIN1(US,RQKL)
ELSE
ARG=-1000.*CWMKL/(US-RQKL)
ARG=AMAX1(ARG,-25.)
CCR(I,LL)= RQKL*(1.-EXP(ARG))
ENDIF
CSMID(I,LL)=AMIN1(US,CCR(I,LL))
!----------------------------------------------------------------------
275 CONTINUE
280 CONTINUE
!----------------------------------------------------------------------
!**********************************************************************
! NOW CHECK THE CLOUDS PRODUCED ABOVE TO MAKE SURE THEY ARE GOOD
! ENOUGH FOR RADIATION CALCULATIONS
!**********************************************************************
!***
!*** NO STRATIFORM CLOUDS FOR THIS TYPE
!***
DO 350 I=MYIS,MYIE
!
LML=LMH(I,J)
LVLIJ=LVL(I,J)
!***
!*** ZERO OUT CLDAMT IF LAND AND BELOW PBOT ABOVE GROUND
!***
! IF(SM(I,J).LT.0.5)THEN
IF((XLAND(I,J)-1.) .LT. 0.5)THEN
DO L=1,LML
LL=LML-L+1+LVLIJ
DDP=PSFC(I)-PMID(I,LL)
IF(DDP.GE.PBOT) GO TO 290
CSMID(I,LL)=0.
ENDDO
290 CONTINUE
ENDIF
!***
!*** CHECK FOR OCEAN STRATUS (LOW CLOUD)
!*** LOOK ONLY OVER OCEAN AND ONLY IF AN INVERSION (DTHDP.LE.-0.05)
!*** IS PRESENT WITH AT LEAST 2 CLOUD FREE LAYERS ABOVE IT
!***
! IF(SM(I,J).GT.0.5)THEN
IF((XLAND(I,J)-1.) .GT. 0.5)THEN
!
!*** FIND BASE OF INVERSION
!
LBASE=LM
DO L=1,LML-1
LL=LML-L+1+LVLIJ
DTHDP=(THMID(I,LL-1)-THMID(I,LL)) &
/(PMID(I,LL-1)-PMID(I,LL))
IF(DTHDP.LE.CLAPSE)THEN
LBASE=LL
GO TO 300
ENDIF
ENDDO
300 CONTINUE
!
!*** CHECK 2 LAYERS ABOVE LBASE FOR DRYNESS
!
IF(CSMID(I,LBASE-1).LE.0..AND.CSMID(I,LBASE-2).LE.0. &
.AND.LBASE.LT.LM)THEN
IF(DTHDP.GT.CLPSE)THEN
CLPFIL=1.-((CLPSE-DTHDP)/DCLPS)
ELSE
CLPFIL=1.
ENDIF
!
DO L=1,LML
LL=LML-L+1+LVLIJ
DDP=PSFC(I)-PMID(I,LL)
IF(DDP.GE.PBOT) GO TO 310
CSMID(I,LL)=CSMID(I,LL)*CLPFIL
ENDDO
310 CONTINUE
!
!*** IF NO INVERSION OR IF CLDS EXIST IN EITHER OF THE 2 LAYERS ABOVE
!*** INVERSION, ZERO OUT CLOUD BELOW PBOT
!
ELSE
DO L=1,LML
LL=LML-L+1+LVLIJ
DDP=PSFC(I)-PMID(I,LL)
IF(DDP.GE.PBOT) GO TO 320
CSMID(I,LL)=0.
ENDDO
320 CONTINUE
!
ENDIF
!------------
ENDIF
!------------
!***
!*** REMOVE HIGH CLOUDS ABOVE THE TROPOPAUSE
!***
L400=LM
DO L=1,LML
LL=LML-L+1+LVLIJ
IF(PMID(I,LL).LE.40000.0)THEN
L400=LL
GO TO 330
ENDIF
ENDDO
330 CONTINUE
!
LTROP=LM
DO LL=L400,2,-1
DTHDP=(THMID(I,LL-1)-THMID(I,LL)) &
/(PMID(I,LL-1)-PMID(I,LL))
!
IF(DTHDP.LT.-0.0025.OR.QMID(I,LL).LE.EPSQ1)THEN
LTROP=LL
GOTO 340
ENDIF
ENDDO
340 IF(LTROP.LT.LM)THEN
DO LL=LTROP,1,-1
CSMID(I,LL)=0.
ENDDO
ENDIF
350 CONTINUE
!
!*********************************************************************
!*****************END OF STRATIFORM CLOUD SECTION*****************
!------------------------CONVECTION--------------------------------
!***
!*** CONVECTIVE CLOUD SECTION
!***
!*** THIS PART WAS MODIFIED TO COMPUTE CONVECTIVE CLOUDS AT EACH
!*** MODEL LAYER BASED ON CONVECTIVE PRECIPITATION RATES. CURRENTLY,
!*** CLOUDS ARE SET TO 0.75*CV(I) BELOW 400MB
!*** AND 0.90*CV(I) ABOVE 400MB TO ACCOUNT FOR CIRRUS CAP
!*** Q.ZHAO 95-3-22
!
!***
!*** NON-PRECIPITATING CLOUD FRACTION OF 20 PERCENT IS ADDED AT
!*** AT POINTS WHERE THE SHALLOW AND DEEP CONVECTIONS ACCUR.
!*** Q. ZHAO 97-5-2
!
! COMPUTE THE CONVECTIVE CLOUD COVER FOR RADIATION
!
!-----------------------------------------------------------------
IF(CNCLD)THEN
!-----------------------------------------------------------------
DO 375 I=MYIS,MYIE
IF(HBOT(I,J)-HTOP(I,J).GT.1.0)THEN
! IF(HTOP(I,J).LT.HBOT(I,J))THEN
SV(I)=0.0
ELSE
SV(I)=0.0
ENDIF
!
PMOD=CUPPT(I,J)*24.0*1000.0/CLSTP
NMOD=0
!
DO NC=1,9
IF(PMOD.GT.PPT(NC)) NMOD=NC
ENDDO
!
!*** CLOUD TOPS AND BOTTOMS COME FROM CUCNVC
!*** ADD LVL TO BE CONSISTENT WITH OTHER WORKING ARRAYS
!
IF(NMOD.EQ.0)THEN
CV(I)=0.
ELSEIF(NMOD.EQ.9)THEN
CV(I)=CC(9)
ELSE
CC1=CC(NMOD)
CC2=CC(NMOD+1)
P1=PPT(NMOD)
P2=PPT(NMOD+1)
CV(I)=CC1+(CC2-CC1)*(PMOD-P1)/(P2-P1)
ENDIF
!
CV(I)=AMAX1(SV(I),CV(I))
CV(I)=AMIN1(1.0,CV(I))
!
IF(CV(I).EQ.0.0)THEN
ICVT(I)=0
ICVB(I)=0
ELSE
ICVT(I)=INT(HTOP(I,J)+0.50)+LVL(I,J)
ICVB(I)=INT(HBOT(I,J)+0.50)+LVL(I,J)
ENDIF
375 CONTINUE
!***
!*** MAKE SURE CLOUDS ARE DEEP ENOUGH
!***
DO I=MYIS,MYIE
BCLD(I)=CV(I).GT.0..AND. (ICVB(I)-ICVT(I)).GE.1
BTEMP1(I)=BCLD(I)
ENDDO
!***
!*** COMPUTE CONVECTIVE CLOUD FRACTION
!***
DO 390 I=MYIS,MYIE
IF(BCLD(I)) THEN
LML=LMH(I,J)
LVLIJ=LVL(I,J)
!
DO L=1,LML
LL=L+LVLIJ
IF(LL.GT.ICVB(I).OR.LL.LT.ICVT(I))THEN
CCMID(I,LL)=0.
ELSE
CCMID(I,LL)=CV(I)
ENDIF
CCMID(I,LL)=AMIN1(1.0,CCMID(I,LL))
ENDDO
ENDIF
390 CONTINUE
!***
!*** REMOVE HIGH CLOUDS ABOVE THE TROPOPAUSE
!***
L400=LM
DO 425 I=MYIS,MYIE
LML=LMH(I,J)
LVLIJ=LVL(I,J)
!
DO L = 1, LML
LL=LML-L+1+LVLIJ
IF(PMID(I,LL).LE.40000.0)THEN
L400=LL
GO TO 400
ENDIF
ENDDO
400 CONTINUE
!
LTROP=LM
DO LL=L400,2,-1
DTHDP=(THMID(I,LL-1)-THMID(I,LL))/(PMID(I,LL-1)-PMID(I,LL))
IF(DTHDP.LT.-0.0025.OR.QMID(I,LL).LE.EPSQ1)THEN
LTROP=LL
GOTO 410
ENDIF
ENDDO
410 IF(LTROP.LT.LM)THEN
DO LL=LTROP,1,-1
CCMID(I,LL)=0.
ENDDO
ENDIF
425 CONTINUE
!***
!-----------------------------------------------------------------
ENDIF
!------------------------CONVECTION--------------------------------
!*****************END OF CONVECTIVE CLOUD SECTION*****************
!*********************************************************************
!***
!*** DETERMINE THE FRACTIONAL CLOUD COVERAGE FOR HIGH, MID
!*** AND LOW OF CLOUDS FROM THE CLOUD COVERAGE AT EACH LEVEL
!***
!*** NOTE: THIS IS FOR DIAGNOSTICS ONLY!!!
!***
!***
! DO 500 I=MYIS,MYIE
!!
! CSTR(I)=0.0
!!
! DO L=0,LM
! CLDAMT(L)=0.
! ENDDO
!!
!!*** NOW GOES LOW, MIDDLE, HIGH
!!
! DO 480 NLVL=1,3
! CLDMAX=0.
! MALVL=LM
! LLTOP=LTOP(NLVL)+LVL(I,J)
!!***
!!*** GO TO THE NEXT CLOUD LAYER IF THE TOP OF THE CLOUD-TYPE IN
!!*** QUESTION IS BELOW GROUND OR IS IN THE LOWEST LAYER ABOVE GROUND.
!!***
! IF(LLTOP.GE.LM)GO TO 480
!!
! IF(NLVL.GT.1)THEN
! LLBOT=LTOP(NLVL-1)-1+LVL(I,J)
! LLBOT=MIN(LLBOT,LM1)
! ELSE
! LLBOT=LM1
! ENDIF
!!
! DO 435 L=LLTOP,LLBOT
! CLDAMT(L)=AMAX1(CSMID(I,L),CCMID(I,L))
! IF(CLDAMT(L).GT.CLDMAX)THEN
! MALVL=L
! CLDMAX=CLDAMT(L)
! ENDIF
! 435 CONTINUE
!!*********************************************************************
!! NOW, CALCULATE THE TOTAL CLOUD FRACTION IN THIS PRESSURE DOMAIN
!! USING THE METHOD DEVELOPED BY Y.H., K.A.C. AND A.K. (NOV., 1992).
!! IN THIS METHOD, IT IS ASSUMED THAT SEPERATED CLOUD LAYERS ARE
!! RADOMLY OVERLAPPED AND ADJACENT CLOUD LAYERS ARE MAXIMUM OVERLAPPED.
!! VERTICAL LOCATION OF EACH TYPE OF CLOUD IS DETERMINED BY THE THICKEST
!! CONTINUING CLOUD LAYERS IN THE DOMAIN.
!!*********************************************************************
! CL1=0.0
! CL2=0.0
! KBT1=LLBOT
! KBT2=LLBOT
! KTH1=0
! KTH2=0
!!
! DO 450 LL=LLTOP,LLBOT
! L=LLBOT-LL+LLTOP
! BIT1=.FALSE.
! CR1=CLDAMT(L)
! BITX=(PINT(I,L).GE.PTOPC(NLVL+1)).AND. &
! (PINT(I,L).LT.PTOPC(NLVL)).AND. &
! (CLDAMT(L).GT.0.0)
! BIT1=BIT1.OR.BITX
! IF(.NOT.BIT1)GO TO 450
!!***
!!*** BITY=T: FIRST CLOUD LAYER; BITZ=T:CONSECUTIVE CLOUD LAYER
!!*** NOTE: WE ASSUME THAT THE THICKNESS OF EACH CLOUD LAYER IN THE
!!*** DOMAIN IS LESS THAN 200 MB TO AVOID TOO MUCH COOLING OR
!!*** HEATING. SO WE SET CTHK(NLVL)=200*E2. BUT THIS LIMIT MAY
!!*** WORK WELL FOR CONVECTIVE CLOUDS. MODIFICATION MAY BE
!!*** NEEDED IN THE FUTURE.
!!***
! BITY=BITX.AND.(KTH2.LE.0)
! BITZ=BITX.AND.(KTH2.GT.0)
!!
! IF(BITY)THEN
! KBT2=L
! KTH2=1
! ENDIF
!!
! IF(BITZ)THEN
! KTOP1=KBT2-KTH2+1
! DPCL=PMID(I,KBT2)-PMID(I,KTOP1)
! IF(DPCL.LT.CTHK(NLVL))THEN
! KTH2=KTH2+1
! ELSE
! KBT2=KBT2-1
! ENDIF
! ENDIF
! IF(BITX)CL2=AMAX1(CL2,CR1)
!!***
!!*** AT THE DOMAIN BOUNDARY OR SEPARATED CLD LAYERS, RANDOM OVERLAP.
!!*** CHOOSE THE THICKEST OR THE LARGEST FRACTION AMT AS THE CLD
!!*** LAYER IN THAT DOMAIN.
!!***
! BIT2=.FALSE.
! BITY=BITX.AND.(CLDAMT(L-1).LE.0.0.OR. &
! PINT(I,L-1).LT.PTOPC(NLVL+1))
! BITZ=BITY.AND.CL1.GT.0.0
! BITW=BITY.AND.CL1.LE.0.0
! BIT2=BIT2.OR.BITY
! IF(.NOT.BIT2)GO TO 450
!!
! IF(BITZ)THEN
! KBT1=INT((CL1*KBT1+CL2*KBT2)/(CL1+CL2))
! KTH1=INT((CL1*KTH1+CL2*KTH2)/(CL1+CL2))+1
! CL1=CL1+CL2-CL1*CL2
! ENDIF
!!
! IF(BITW)THEN
! KBT1=KBT2
! KTH1=KTH2
! CL1=CL2
! ENDIF
!!
! IF(BITY)THEN
! KBT2=LLBOT
! KTH2=0
! CL2=0.0
! ENDIF
! 450 CONTINUE
!!***
! CLDCFR(I,NLVL)=AMIN1(1.0,CL1)
! MTOP(I,NLVL)=MIN(KBT1,KBT1-KTH1+1)
! MBOT(I,NLVL)=KBT1
! 480 CONTINUE
! 500 CONTINUE
!***
!*** SET THE UN-NEEDED TAUDAR TO ONE
!***
DO I=MYIS,MYIE
TAUDAR(I)=1.0
ENDDO
!----------------------------------------------------------------------
! NOW, CALCULATE THE CLOUD RADIATIVE PROPERTIES AFTER DAVIS (1982),
! HARSHVARDHAN ET AL (1987) AND Y.H., K.A.C. AND A.K. (1993).
!
! UPDATE: THE FOLLOWING PARTS ARE MODIFIED, AFTER Y.T.H. (1994), TO
! CALCULATE THE RADIATIVE PROPERTIES OF CLOUDS ON EACH MODEL
! LAYER. BOTH CONVECTIVE AND STRATIFORM CLOUDS ARE USED
! IN THIS CALCULATIONS.
!
! QINGYUN ZHAO 95-3-22
!
!----------------------------------------------------------------------
!
!***
!*** INITIALIZE ARRAYS FOR USES LATER
!***
DO 600 I=MYIS,MYIE
LML=LMH(I,J)
LVLIJ=LVL(I,J)
!
!***
!*** NOTE: LAYER=1 IS THE SURFACE, AND LAYER=2 IS THE FIRST CLOUD
!*** LAYER ABOVE THE SURFACE AND SO ON.
!***
EMIS(I,1)=1.0
KTOP(I,1)=LP1
KBTM(I,1)=LP1
CAMT(I,1)=1.0
ITYP(I,1)=0
KCLD(I)=2
!
DO NBAND=1,NB
RRCL(I,NBAND,1)=0.0
TTCL(I,NBAND,1)=1.0
ENDDO
!
DO 510 L=2,LP1
ITYP(I,L)=0
CAMT(I,L)=0.0
KTOP(I,L)=1
KBTM(I,L)=1
EMIS(I,L)=0.0
!
DO NBAND=1,NB
RRCL(I,NBAND,L)=0.0
TTCL(I,NBAND,L)=1.0
ENDDO
510 CONTINUE
!***
!*** NOW CALCULATE THE AMOUNT, TOP, BOTTOM AND TYPE OF EACH CLOUD LAYER
!*** CLOUD TYPE=1: STRATIFORM CLOUD
!*** TYPE=2: CONVECTIVE CLOUD
!*** WHEN BOTH CONVECTIVE AND STRATIFORM CLOUDS EXIST AT THE SAME POINT,
!*** SELECT CONVECTIVE CLOUD (TYPE=2),IN OTHER WORDS, CONVECTIVE CLOUDS
!*** HAVE THE HIGHER PRIORITY THAN STRATIFORM CLOUDS.
!*** CLOUD LAYERS ARE SEPARATED BY:
!*** 1. NO-CLOUD LAYER
!*** 2. DIFFERENT CLOUD TYPE
!*** NOTE: THERE IS ONLY ONE CONVECTIVE CLOUD LAYER IN ONE COLUMN.
!*** KTOP AND KBTM ARE THE TOP AND BOTTOM OF EACH CLOUD LAYER IN TERMS O
!*** ETA MODEL LEVEL.
!***
DO 540 L=2,LML
LL=LML-L+1+LVLIJ
BITC=CCMID(I,LL).GT.0.1
BITS=CSMID(I,LL).GT.0.1
BITCP1=CCMID(I,LL+1).GT.0.1
BITSP1=CSMID(I,LL+1).GT.0.1
BIT1=BITS.OR.BITC
!-------------------
IF(BIT1)THEN
!-------------------
IF(ITYP(I,KCLD(I)).EQ.0)THEN
CAMT(I,KCLD(I))=CSMID(I,LL)
ITYP(I,KCLD(I))=1
KBTM(I,KCLD(I))=LL
!
IF(BITC)THEN
CAMT(I,KCLD(I))=CCMID(I,LL)
ITYP(I,KCLD(I))=2
ENDIF
ELSE
IF(BITC)THEN
IF(BITCP1)THEN
CAMT(I,KCLD(I))=AMAX1(CAMT(I,KCLD(I)),CCMID(I,LL))
ELSE
KCLD(I)=KCLD(I)+1
CAMT(I,KCLD(I))=CCMID(I,LL)
ITYP(I,KCLD(I))=2
KTOP(I,KCLD(I)-1)=LL+1
KBTM(I,KCLD(I))=LL
ENDIF
ELSE
IF(BITCP1)THEN
KCLD(I)=KCLD(I)+1
CAMT(I,KCLD(I))=CSMID(I,LL)
ITYP(I,KCLD(I))=1
KTOP(I,KCLD(I)-1)=LL+1
KBTM(I,KCLD(I))=LL
ELSE
CAMT(I,KCLD(I))=AMAX1(CAMT(I,KCLD(I)),CSMID(I,LL))
ENDIF
ENDIF
ENDIF
!-------------------
ELSE
!-------------------
IF(BITCP1.OR.BITSP1)THEN
KCLD(I)=KCLD(I)+1
KTOP(I,KCLD(I)-1)=LL+1
ITYP(I,KCLD(I))=0
CAMT(I,KCLD(I))=0.0
ENDIF
!-------------------
ENDIF
!-------------------
540 CONTINUE
!***
!*** THE REAL NUMBER OF CLOUD LAYERS IS (THE FIRST IS THE GROUNG;
!*** THE LAST IS THE SKY):
!***
NCLDS(I)=KCLD(I)-2
NCLD=NCLDS(I)
!***
!*** NOW CALCULATE CLOUD RADIATIVE PROPERTIES
!***
IF(NCLD.GE.1)THEN
!***
!*** NOTE: THE FOLLOWING CALCULATIONS, THE UNIT FOR PRESSURE IS MB!!!
!***
DO 580 NC=2,NCLD+1
!
TAUC(I)=0.0
QSUM=0.0
NKTP=LP1
NBTM=0
BITX=CAMT(I,NC).GT.0.1
NKTP=MIN(NKTP,KTOP(I,NC))
NBTM=MAX(NBTM,KBTM(I,NC))
!
DO 560 LL=NKTP,NBTM
IF(LL.GE.KTOP(I,NC).AND.LL.LE.KBTM(I,NC).AND.BITX)THEN
PRS1=PINT(I,LL)*0.01
PRS2=PINT(I,LL+1)*0.01
DELP=PRS2-PRS1
TCLD=TMID(I,LL)-273.16
QSUM=QSUM+QMID(I,LL)*DELP*(PRS1+PRS2) &
/(120.1612*SQRT(TMID(I,LL)))
!***
!*** FOR CONVECTIVE CLOUD OR STARTIFORM CLOUD WITH TOP ABOVE 500MB
!***
IF(ITYP(I,NC).EQ.2 &
.OR.PINT(I,KTOP(I,NC)).LE.PTOPC(3))THEN
IF(TCLD.LE.-10.0)THEN
TAUC(I)=TAUC(I)+DELP*AMAX1(0.1E-3, &
2.0E-6*(TCLD+82.5)**2)
ELSE
TAUC(I)=TAUC(I)+DELP*AMIN1(0.08,6.949E-3*TCLD+0.1)
ENDIF
ELSE
!***
!*** FOR LOW AND MID STRATIFORM CLOUDS
!***
IF(TCLD.LE.-20.0)THEN
TAUC(I)=TAUC(I)+DELP*AMAX1(0.1E-3,2.56E-5* &
(TCLD+82.5)**2)
ELSE
TAUC(I)=TAUC(I)+DELP*0.1
ENDIF
ENDIF
ENDIF
560 CONTINUE
!
IF(BITX)EMIS(I,NC)=1.0-EXP(-0.75*TAUC(I))
IF(QSUM.GE.EPSQ1)THEN
!
DO 570 NBAND=1,NB
IF(BITX)THEN
PROD=ABCFF(NBAND)*QSUM
DDX=TAUC(I)/(TAUC(I)+PROD)
EEX=1.0-DDX
IF(ABS(EEX).GE.1.E-8)THEN
DD=DDX
EE=EEX
FF=1.0-DD*0.85
AA=MIN(50.0,SQRT(3.0*EE*FF)*TAUC(I))
AA=EXP(-AA)
BB=FF/EE
GG=SQRT(BB)
DD=(GG+1.0)*(GG+1.0)-(GG-1.0)*(GG-1.0)*AA*AA
RRCL(I,NBAND,NC)=MAX(0.1E-5,(BB-1.0)*(1.0-AA*AA)/DD)
TTCL(I,NBAND,NC)=AMAX1(0.1E-5,4.0*GG*AA/DD)
ENDIF
ENDIF
570 CONTINUE
ENDIF
580 CONTINUE
!
ENDIF
!
600 CONTINUE
!*********************************************************************
!****************** COMPUTE OZONE AT MIDLAYERS *********************
!*********************************************************************
!
!*** MODIFY PRESSURES SO THAT THE ENTIRE COLUMN OF OZONE (TO 0 MB)
!*** IS INCLUDED IN THE MODEL COLUMN EVEN WHEN PT > 0 MB
!***
DO L=1,LM
DO I=MYIS,MYIE
DENOM=1./(PINT(I,LP1)-PINT(I,1))
FCTRA=PINT(I,LP1)*DENOM
FCTRB=-PINT(I,1)*PINT(I,LP1)*DENOM
POZN(I,L)=PMID(I,L)*FCTRA+FCTRB
ENDDO
ENDDO
!
! CALL OZON2D(LM,POZN,XLAT,RSIN1,RCOS1,RCOS2,OZN)
CALL OZON2D(LM,POZN,XLAT,RSIN1,RCOS1,RCOS2,OZN, &
! XDUO3N,XDO3N4,XDO3N2,XDO3N3, &
! PRGFDL,MYIS,MYIE, &
MYIS,MYIE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!
!***
!*** NOW THE VARIABLES REQUIRED BY RADFS HAVE BEEN CALCULATED.
!***
!----------------------------------------------------------------------
!***
!*** CALL THE GFDL RADIATION DRIVER
!***
!***
CALL RADFS &
(PSFC,PMID,PINT,QMID,TMID,OZN,TSKN,SLMSK,ALBEDO,XLAT &
, CAMT,KTOP,KBTM,NCLDS,EMIS,RRCL,TTCL &
, COSZ,TAUDAR,1 &
, 1,0 &
! , ETAD,AETA,ITIMSW,ITIMLW,JD,R1,HOUR,TENDS,TENDL &
, ITIMSW,ITIMLW,JD,HOUR &
, TENDS(its,kts,j),TENDL(its,kts,j) &
! , T1,T2,T4,EM1V,EM1VW,EM3V &
! , TABLE1,TABLE2,TABLE3,EM1,EM1WDE,EM3 &
, FLWUP,FSWUP,FSWDN,FSWDNS,FSWUPS,FLWDNS,FLWUPS &
! , RAD1,RAD2,RAD3,RAD4 &
! , SKO3R,AB15WD,SKC1R,SKO2D &
, ids,ide, jds,jde, kds,kde &
, ims,ime, jms,jme, kms,kme &
, its,ite, jts,jte, kts,kte )
!----------------------------------------------------------------------
DO I=MYIS,MYIE
GLW(I,J) = FLWDNS(I)
GSW(I,J) = FSWDNS(I)-FSWUPS(I)
ENDDO
DO 650 I=MYIS,MYIE
! PDSLIJ=PDSL(I,J)
PMOD=CUPPT(I,J)*24.0*1000.0/CLSTP
! CFRACL(I,J)=CLDCFR(I,1)
! CFRACM(I,J)=CLDCFR(I,2)
! CFRACH(I,J)=CLDCFR(I,3)
!
!*** ARRAYS ACFRST AND ACFRCV ACCUMULATE AVERAGE STRATIFORM AND
!*** CONVECTIVE CLOUD FRACTIONS, RESPECTIVELY. THIS INFORMATION
!*** IS PASSED TO THE POST PROCESSOR VIA COMMON BLOCK ACMCLD.
!
!SH CFRAVG=AMAX1(CFRACL(I,J),AMAX1(CFRACM(I,J),CFRACH(I,J)))
IF(CNCLD)THEN
IF(PMOD.LE.PPT(1))THEN
!SH ACFRST(I,J)=ACFRST(I,J)+CFRAVG
!SH NCFRST(I,J)=NCFRST(I,J)+1
ELSE
!SH ACFRCV(I,J)=ACFRCV(I,J)+CFRAVG
!SH NCFRCV(I,J)=NCFRCV(I,J)+1
ENDIF
ELSE
!SH ACFRST(I,J)=ACFRST(I,J)+CFRAVG
!SH NCFRST(I,J)=NCFRST(I,J)+1
ENDIF
650 CONTINUE
!***
!*** COLLECT ATMOSPHERIC TEMPERATURE TENDENCIES DUE TO RADIATION.
!*** ALSO COLLECT THE TOTAL SW AND INCOMING LW RADIATION (W/M**2)
!*** AND CONVERT TO FORM NEEDED FOR PREDICTION OF THS IN SURFCE.
!***
DO 660 I=MYIS,MYIE
DO L=1,LM
LL=LVL(I,J)+L
IF(SHORT)RSWTT(I,L,J)=TENDS(I,LL,J)
IF(LONG) RLWTT(I,L,J)=TENDL(I,LL,J)
IF(LL.EQ.LM)GO TO 660
ENDDO
660 CONTINUE
!***
!*** SUM THE LW INCOMING AND SW RADIATION (W/M**2) FOR RADIN.
!***
DO 675 I=MYIS,MYIE
IF(LONG)THEN
SIGT4(I,J)=STBOL*TMID(I,LM)*TMID(I,LM)* &
TMID(I,LM)*TMID(I,LM)
ENDIF
!
!*** ACCUMULATE VARIOUS LW AND SW RADIATIVE FLUXES FOR POST
!*** PROCESSOR. PASSED VIA COMMON ACMRDL AND ACMRDS.
!
IF(LONG)THEN
RLWIN(I,J) =FLWDNS(I)
RLWOUT(I,J)=FLWUPS(I)
RLWTOA(I,J)=FLWUP(I)
ENDIF
IF(SHORT)THEN
RSWIN(I,J) =FSWDNS(I)
RSWOUT(I,J)=FSWUPS(I)
RSWTOA(I,J)=FSWUP(I)
ENDIF
675 CONTINUE
!***
!*** THIS ROW IS FINISHED. GO TO NEXT
!***
! *********************
700 CONTINUE
! *********************
!----------------------------------------------------------------------
!***
!*** CALLS TO RADIATION THIS TIME STEP ARE COMPLETE.
!***
!----------------------------------------------------------------------
!----------------------------------------------------------------------
!***
!*** HORIZONTAL SMOOTHING OF TEMPERATURE TENDENCIES
!***
!----------------------------------------------------------------------
IF(SHORT) THEN
DO 800 L=1,LM
CALL ZERO2(TL, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
CALL ZERO2(FNE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
CALL ZERO2(FSE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!
! IF(KSMUD.GE.1)THEN
! DO 750 KS=1,KSMUD
!!
! DO J=MYJS,MYJE
! DO I=MYIS,MYIE
! TL(I,J)=RSWTT(I,L,J)
!! TL(I,J)=RSWTT(I,L,J)*HTM(I,L,J)
! ENDDO
! ENDDO
!!
! DO J=MYJS,MYJE
! DO I=MYIS,MYIE
! FNE(I,J)=(TL(I+IHE(J),J+1)-TL(I,J))
!! *HTM(I,L,J)*HTM(I+IHE(J),J+1,L)
! ENDDO
! ENDDO
!!
! DO J=MYJS1,MYJE
! DO I=MYIS,MYIE
! FSE(I,J)=(TL(I+IHE(J),J-1)-TL(I,J))
!! *HTM(I+IHE(J),J-1,L)*HTM(I,L,J)
! ENDDO
! ENDDO
!!
! DO J=MYJS2,MYJE2
! DO I=MYIS,MYIE
! TL(I,J)=(FNE(I,J)-FNE(I+IHW(J),J-1) &
! +FSE(I,J)-FSE(I+IHW(J),J+1)) &
! *0.125+TL(I,J)
!! *HBM2(I,J)*0.125+TL(I,J)
! ENDDO
! ENDDO
!!
! DO J=MYJS,MYJE
! DO I=MYIS,MYIE
! RSWTT(I,L,J)=TL(I,J)
! ENDDO
! ENDDO
!!
! 750 CONTINUE
! ENDIF
!
800 CONTINUE
ENDIF
!----------------------------------------------------------------------
!
IF(LONG)THEN
!
DO 900 L=1,LM
CALL ZERO2(TL, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
CALL ZERO2(FNE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
CALL ZERO2(FSE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!
! IF(KSMUD.GE.1)THEN
! DO 850 KS=1,KSMUD
!!
! DO J=MYJS,MYJE
! DO I=MYIS,MYIE
! TL(I,J)=RLWTT(I,L,J)
!! TL(I,J)=RLWTT(I,L,J)*HTM(I,L,J)
! ENDDO
! ENDDO
!!
! DO J=MYJS,MYJE1
! DO I=MYIS,MYIE
! FNE(I,J)=(TL(I+IHE(J),J+1)-TL(I,J))
!! *HTM(I,L,J)*HTM(I+IHE(J),J+1,L)
! ENDDO
! ENDDO
!!
! DO J=MYJS1,MYJE
! DO I=MYIS,MYIE
! FSE(I,J)=(TL(I+IHE(J),J-1)-TL(I,J))
!! *HTM(I+IHE(J),J-1,L)*HTM(I,L,J)
! ENDDO
! ENDDO
!!
! DO J=MYJS2,MYJE2
! DO I=MYIS,MYIE
! TL(I,J)=(FNE(I,J)-FNE(I+IHW(J),J-1) &
! +FSE(I,J)-FSE(I+IHW(J),J+1)) &
! *0.125+TL(I,J)
!! *HBM2(I,J)*0.125+TL(I,J)
! ENDDO
! ENDDO
!!
! DO J=MYJS,MYJE
! DO I=MYIS,MYIE
! RLWTT(I,L,J)=TL(I,J)
! ENDDO
! ENDDO
!!
! 850 CONTINUE
! ENDIF
900 CONTINUE
ENDIF
!-----------------------------------------------------------------------
!***
!*** RESET CUPPT,HTOP,HBOT FROM THIS CALL TO RADTN
!***
IF(MOD(NTSD,NRADPP).EQ.1)THEN
DO J=MYJS,MYJE
DO I=MYIS,MYIE
CUPPT(I,J)=0.
HTOP(I,J)=100.
HBOT(I,J)=0.
ENDDO
ENDDO
ENDIF
!-----------------------------------------------------------------------
END SUBROUTINE RADTN
!-----------------------------------------------------------------------
SUBROUTINE ZENITH(TIMES,DAYI,HOUR,IDAT,IHRST,GLON,GLAT,CZEN, & 2
MYIS,MYIE,MYJS,MYJE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!----------------------------------------------------------------------
IMPLICIT NONE
!----------------------------------------------------------------------
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
INTEGER, INTENT(IN) :: MYJS,MYJE,MYIS,MYIE
REAL, INTENT(IN) :: TIMES
REAL, INTENT(OUT) :: HOUR,DAYI
INTEGER, INTENT(IN) :: IHRST
INTEGER, INTENT(IN), DIMENSION(3) :: IDAT
REAL, INTENT(IN ), DIMENSION(its:ite,jts:jte) :: GLAT,GLON
REAL, INTENT(OUT), DIMENSION(its:ite,jts:jte) :: CZEN
REAL, PARAMETER :: GSTC1=24110.54841,GSTC2=8640184.812866, &
GSTC3=9.3104E-2,GSTC4=-6.2E-6, &
PI=3.1415926,PI2=2.*PI,PIH=0.5*PI, &
DEG2RD=1.745329E-2,OBLIQ=23.440*DEG2RD, &
ZEROJD=2451545.0
REAL :: DAY,YFCTR,ADDDAY,STARTYR,DATJUL,DIFJD,SLONM, &
ANOM,SLON,DEC,RA,DATJ0,TU,STIM0,SIDTIM,HRANG
REAL :: HRLCL,SINALT
INTEGER :: KMNTH,KNT,IDIFYR,J,I
LOGICAL :: LEAP
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
INTEGER :: MONTH (12)
!-----------------------------------------------------------------------
DATA MONTH/31,28,31,30,31,30,31,31,30,31,30,31/
!***********************************************************************
! SAVE MONTH
DAY=0.
LEAP=.FALSE.
IF(MOD(IDAT(3),4).EQ.0)THEN
MONTH(2)=29
LEAP=.TRUE.
ENDIF
IF(IDAT(1).GT.1)THEN
KMNTH=IDAT(1)-1
DO 10 KNT=1,KMNTH
DAY=DAY+REAL(MONTH(KNT))
10 CONTINUE
ENDIF
!***
!*** CALCULATE EXACT NUMBER OF DAYS FROM BEGINNING OF YEAR TO
!*** FORECAST TIME OF INTEREST
!***
DAY=DAY+REAL(IDAT(2)-1)+(REAL(IHRST)+TIMES/3600.)/24.
DAYI=REAL(INT(DAY)+1)
HOUR=(DAY-DAYI+1.)*24.
YFCTR=2000.-IDAT(3)
!-----------------------------------------------------------------------
!***
!*** FIND CELESTIAL LONGITUDE OF THE SUN THEN THE SOLAR DECLINATION AND
!*** RIGHT ASCENSION.
!***
!-----------------------------------------------------------------------
IDIFYR=IDAT(3)-2000
!***
!*** FIND JULIAN DATE OF START OF THE RELEVANT YEAR
!*** ADDING IN LEAP DAYS AS NEEDED
!***
IF(IDIFYR.LT.0)THEN
ADDDAY=REAL(IDIFYR/4)
ELSE
ADDDAY=REAL((IDIFYR+3)/4)
ENDIF
STARTYR=ZEROJD+IDIFYR*365.+ADDDAY-0.5
!***
!*** THE JULIAN DATE OF THE TIME IN QUESTION
!***
DATJUL=STARTYR+DAY
!
!*** DIFFERENCE OF ACTUAL JULIAN DATE FROM JULIAN DATE
!*** AT 00H 1 January 2000
!
DIFJD=DATJUL-ZEROJD
!
!*** MEAN GEOMETRIC LONGITUDE OF THE SUN
!
SLONM=(280.460+0.9856474*DIFJD)*DEG2RD+YFCTR*PI2
!
!*** THE MEAN ANOMOLY
!
ANOM=(357.528+0.9856003*DIFJD)*DEG2RD
!
!*** APPARENT GEOMETRIC LONGITUDE OF THE SUN
!
SLON=SLONM+(1.915*SIN(ANOM)+0.020*SIN(2.*ANOM))*DEG2RD
IF(SLON.GT.PI2)SLON=SLON-PI2
!
!*** DECLINATION AND RIGHT ASCENSION
!
DEC=ASIN(SIN(SLON)*SIN(OBLIQ))
RA=ACOS(COS(SLON)/COS(DEC))
IF(SLON.GT.PI)RA=PI2-RA
!***
!*** FIND THE GREENWICH SIDEREAL TIME THEN THE LOCAL SOLAR
!*** HOUR ANGLE.
!***
DATJ0=STARTYR+DAYI-1.
TU=(DATJ0-2451545.)/36525.
STIM0=GSTC1+GSTC2*TU+GSTC3*TU**2+GSTC4*TU**3
SIDTIM=STIM0/3600.+YFCTR*24.+1.00273791*HOUR
SIDTIM=SIDTIM*15.*DEG2RD
IF(SIDTIM.LT.0.)SIDTIM=SIDTIM+PI2
IF(SIDTIM.GT.PI2)SIDTIM=SIDTIM-PI2
HRANG=SIDTIM-RA
!
DO 100 J=MYJS,MYJE
DO 100 I=MYIS,MYIE
! HRLCL=HRANG-GLON(I,J)
HRLCL=HRANG+GLON(I,J)+PI2
!***
!*** THE ZENITH ANGLE IS THE COMPLEMENT OF THE ALTITUDE THUS THE
!*** COSINE OF THE ZENITH ANGLE EQUALS THE SINE OF THE ALTITUDE.
!***
SINALT=SIN(DEC)*SIN(GLAT(I,J))+COS(DEC)*COS(HRLCL)* &
COS(GLAT(I,J))
IF(SINALT.LT.0.)SINALT=0.
CZEN(I,J)=SINALT
100 CONTINUE
!***
!*** IF THE FORECAST IS IN A DIFFERENT YEAR THAN THE START TIME,
!*** RESET DAYI TO THE PROPER DAY OF THE NEW YEAR (IT MUST NOT BE
!*** RESET BEFORE THE SOLAR ZENITH ANGLE IS COMPUTED).
!***
IF(DAYI.GT.365.)THEN
IF(.NOT.LEAP)THEN
DAYI=DAYI-365.
ELSEIF(LEAP.AND.DAYI.GT.366.)THEN
DAYI=DAYI-366.
ENDIF
ENDIF
END SUBROUTINE ZENITH
!-----------------------------------------------------------------------
SUBROUTINE OZON2D (LK,POZN,XLAT,RSIN1,RCOS1,RCOS2,QO3, & 1
! XDUO3N,XDO3N4,XDO3N2,XDO3N3, &
! PRGFDL,MYIS,MYIE, &
MYIS,MYIE, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!----------------------------------------------------------------------
IMPLICIT NONE
!----------------------------------------------------------------------
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
INTEGER, INTENT(IN) :: LK,MYIS,MYIE
REAL, INTENT(IN), DIMENSION(its:ite,kts:kte) :: POZN
REAL, INTENT(IN), DIMENSION(its:ite) :: XLAT
REAL, INTENT(INOUT), DIMENSION(its:ite,kts:kte) :: QO3
REAL, INTENT(IN) :: RSIN1,RCOS1,RCOS2
!----------------------------------------------------------------------
INTEGER, PARAMETER :: NL=81,NLP1=NL+1,LNGTH=37*NL
REAL, PARAMETER :: RTD=57.2957795
! REAL, INTENT(IN), DIMENSION(37,NL) :: XDUO3N,XDO3N4,XDO3N2,XDO3N3
! REAL, INTENT(IN), DIMENSION(NL) :: PRGFDL
!----------------------------------------------------------------------
!----------------------------------------------------------------------
INTEGER,DIMENSION(its:ite) :: JJROW
REAL, DIMENSION(its:ite) :: TTHAN
REAL, DIMENSION(its:ite,NL) :: QO3O3
INTEGER :: I,K,NUMITR,ILOG,IT,NHALF
REAL :: TH2,DO3V,DO3VP,APHI,APLO
!----------------------------------------------------------------------
DO I=MYIS,MYIE
TH2=0.2*XLAT(I)
JJROW(I)=19.001-TH2
TTHAN(I)=(19-JJROW(I))-TH2
ENDDO
!
!*** SEASONAL AND SPATIAL INTERPOLATION DONE BELOW.
!
DO K=1,NL
DO I=MYIS,MYIE
DO3V=XDUO3N(JJROW(I),K)+RSIN1*XDO3N2(JJROW(I),K) &
+RCOS1*XDO3N3(JJROW(I),K) &
+RCOS2*XDO3N4(JJROW(I),K)
DO3VP=XDUO3N(JJROW(I)+1,K)+RSIN1*XDO3N2(JJROW(I)+1,K) &
+RCOS1*XDO3N3(JJROW(I)+1,K) &
+RCOS2*XDO3N4(JJROW(I)+1,K)
!
!*** NOW LATITUDINAL INTERPOLATION
!*** AND CONVERT O3 INTO MASS MIXING RATIO (ORIG DATA MPY BY 1.E4)
!
QO3O3(I,K)=1.E-4*(DO3V+TTHAN(I)*(DO3VP-DO3V))
ENDDO
ENDDO
!***
!*** VERTICAL INTERPOLATION FOR EACH GRIDPOINT (LINEAR IN LN P)
!***
NUMITR=0
ILOG=NL
20 CONTINUE
ILOG=(ILOG+1)/2
IF(ILOG.EQ.1)GO TO 25
NUMITR=NUMITR+1
GO TO 20
25 CONTINUE
!
DO 60 K=1,LK
!
NHALF=(NL+1)/2
DO I=MYIS,MYIE
JJROW(I)=NHALF
ENDDO
!
DO 40 IT=1,NUMITR
NHALF=(NHALF+1)/2
DO I=MYIS,MYIE
IF(POZN(I,K).LT.PRGFDL(JJROW(I)-1))THEN
JJROW(I)=JJROW(I)-NHALF
ELSEIF(POZN(I,K).GE.PRGFDL(JJROW(I)))THEN
JJROW(I)=JJROW(I)+NHALF
ENDIF
JJROW(I)=MIN(JJROW(I),NL)
JJROW(I)=MAX(JJROW(I),2)
ENDDO
40 CONTINUE
!
DO 50 I=MYIS,MYIE
IF(POZN(I,K).LT.PRGFDL(1))THEN
QO3(I,K)=QO3O3(I,1)
ELSE IF(POZN(I,K).GT.PRGFDL(NL))THEN
QO3(I,K)=QO3O3(I,NL)
ELSE
APLO=ALOG(PRGFDL(JJROW(I)-1))
APHI=ALOG(PRGFDL(JJROW(I)))
QO3(I,K)=QO3O3(I,JJROW(I))+(ALOG(POZN(I,K))-APHI)/ &
(APLO-APHI)* &
(QO3O3(I,JJROW(I)-1)-QO3O3(I,JJROW(I)))
ENDIF
50 CONTINUE
!
60 CONTINUE
END SUBROUTINE OZON2D
!-----------------------------------------------------------------------
SUBROUTINE ZERO2(ARRAY, & 6
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!----------------------------------------------------------------------
IMPLICIT NONE
!----------------------------------------------------------------------
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
REAL, INTENT(INOUT), DIMENSION(its:ite,jts:jte) :: ARRAY
INTEGER :: I,J
!----------------------------------------------------------------------
DO J=jts,jte
DO I=its,ite
ARRAY(I,J)=0.
ENDDO
ENDDO
END SUBROUTINE ZERO2
!----------------------------------------------------------------
SUBROUTINE O3INT(PHALF,DDUO3N,DDO3N2,DDO3N3,DDO3N4, & 1
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!----------------------------------------------------------------------
IMPLICIT NONE
!----------------------------------------------------------------------
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde , &
ims,ime, jms,jme, kms,kme , &
its,ite, jts,jte, kts,kte
!FPP$ NOCONCUR R
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
! . . . .
! SUBPROGRAM: O3INT COMPUTE ZONAL MEAN OZONE FOR ETA LYRS
! PRGMMR: KENNETH CAMPANA ORG: W/NMC23 DATE: 89-07-07
! MICHAEL BALDWIN ORG: W/NMC22 DATE: 92-06-08
!
! ABSTRACT: THIS CODE WRITTEN AT GFDL...
! CALCULATES SEASONAL ZONAL MEAN OZONE,EVERY 5 DEG OF LATITUDE,
! FOR CURRENT MODEL VERTICAL COORDINATE. OUTPUT DATA IN G/G * 1.E4
! CODE IS CALLED ONLY ONCE.
!
! PROGRAM HISTORY LOG:
! 84-01-01 FELS AND SCHWARZKOPF,GFDL.
! 89-07-07 K. CAMPANA - ADAPTED STAND-ALONE CODE FOR IN-LINE USE.
! 92-06-08 M. BALDWIN - UPDATE TO RUN IN ETA MODEL
!
! USAGE: CALL O3INT(O3,SIGL) OLD
! INPUT ARGUMENT LIST:
! PHALF - MID LAYER PRESSURE (K=LM+1 IS MODEL SURFACE)
! OUTPUT ARGUMENT LIST:
! DDUO3N - ZONAL MEAN OZONE DATA IN ALL MODEL LAYERS (G/G*1.E4)
! DDO3N2 DIMENSIONED(L,N),WHERE L(=37) IS LATITUDE BETWEEN
! DDO3N3 N AND S POLES,N=NUM OF VERTICAL LYRS(K=1 IS TOP LYR)
! DDO3N4 AND SEASON-WIN,SPR,SUM,FALL.
! IN COMMON
!
! OUTPUT FILES:
! OUTPUT - PRINT FILE.
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 200.
!
!$$$
!.... PROGRAM O3INT FROM DAN SCHWARZKOPF-GETS ZONAL MEAN O3
!.. OUTPUT O3 IS WINTER,SPRING,SUMMER,FALL (NORTHERN HEMISPHERE)
!-----------------------------------------------------------------------
! INCLUDE "parmeta"
!-----------------------------------------------------------------------
! *********************************************************
INTEGER :: N,NP,NP2,NM1
! PARAMETER (N=LM,NP=N+1,NP2=N+2,NM1=N-1)
! *********************************************************
!-----------------------------------------------------------------------
!***
!*** SEASONAL CLIMATOLOGIES OF O3 (OBTAINED FROM A PREVIOUSLY RUN
!*** CODE WHICH INTERPOLATES O3 TO USER VERTICAL COORDINATE).
!*** DEFINED AS 5 DEG LAT MEANS N.P.->S.P.
!***
REAL, INTENT(OUT), DIMENSION(37,kte):: DDUO3N,DDO3N2,DDO3N3,DDO3N4
! C O M M O N /SAVMEM/
! ...WINTER.... ...SPRING.... ...SUMMER.... ....FALL.....
! 1 DDUO3N(37,LM), DDO3N2(37,LM), DDO3N3(37,LM), DDO3N4(37,LM)
! ..... K.CAMPANA OCTOBER 1988
!CCC DIMENSION T41(NP2,2),O3O3(37,N,4)
! DIMENSION SIGL(N)
! *********************************************************
REAL :: QI(82)
REAL :: DDUO3(19,kts:kte),RO31(10,41),RO32(10,41),DUO3N(19,41)
REAL :: TEMPN(19)
REAL :: O3HI(10,25),O3LO1(10,16),O3LO2(10,16),O3LO3(10,16), &
O3LO4(10,16)
REAL :: O3HI1(10,16),O3HI2(10,9),PH1(45),PH2(37),P1(48),P2(33)
REAL :: O35DEG(37,kts:kte)
REAL :: RSTD(81),RO3(10,41),RO3M(10,40),RBAR(kts:kte),RDATA(81), &
PHALF(kts:kte+1),P(81),PH(82)
INTEGER :: NKK,NK,NKP,K,L,NCASE,ITAPE,IPLACE,NKMM,NKM,KI,KK,KQ,JJ,KEN
REAL :: O3RD,O3TOT,O3DU
! EQUIVALENCE (O3HI1(1,1),O3HI(1,1)),(O3HI2(1,1),O3HI(1,17))
! EQUIVALENCE (PH1(1),PH(1)),(PH2(1),PH(46))
! EQUIVALENCE (P1(1),P(1)),(P2(1),P(49))
DATA PH1/ 0., &
0.1027246E-04, 0.1239831E-04, 0.1491845E-04, 0.1788053E-04, &
0.2135032E-04, 0.2540162E-04, 0.3011718E-04, 0.3558949E-04, &
0.4192172E-04, 0.4922875E-04, 0.5763817E-04, 0.6729146E-04, &
0.7834518E-04, 0.9097232E-04, 0.1053635E-03, 0.1217288E-03, &
0.1402989E-03, 0.1613270E-03, 0.1850904E-03, 0.2119495E-03, &
0.2423836E-03, 0.2768980E-03, 0.3160017E-03, 0.3602623E-03, &
0.4103126E-03, 0.4668569E-03, 0.5306792E-03, 0.6026516E-03, &
0.6839018E-03, 0.7759249E-03, 0.8803303E-03, 0.9987843E-03, &
0.1133178E-02, 0.1285955E-02, 0.1460360E-02, 0.1660001E-02, &
0.1888764E-02, 0.2151165E-02, 0.2452466E-02, 0.2798806E-02, &
0.3197345E-02, 0.3656456E-02, 0.4185934E-02, 0.4797257E-02/
DATA PH2/ &
0.5503893E-02, 0.6321654E-02, 0.7269144E-02, 0.8368272E-02, &
0.9644873E-02, 0.1112946E-01, 0.1285810E-01, 0.1487354E-01, &
0.1722643E-01, 0.1997696E-01, 0.2319670E-01, 0.2697093E-01, &
0.3140135E-01, 0.3660952E-01, 0.4274090E-01, 0.4996992E-01, &
0.5848471E-01, 0.6847525E-01, 0.8017242E-01, 0.9386772E-01, &
0.1099026E+00, 0.1286765E+00, 0.1506574E+00, 0.1763932E+00, &
0.2065253E+00, 0.2415209E+00, 0.2814823E+00, 0.3266369E+00, &
0.3774861E+00, 0.4345638E+00, 0.4984375E+00, 0.5697097E+00, &
0.6490189E+00, 0.7370409E+00, 0.8344896E+00, 0.9421190E+00, &
0.1000000E+01/
DATA P1/ &
0.9300000E-05, 0.1129521E-04, 0.1360915E-04, 0.1635370E-04, &
0.1954990E-04, 0.2331653E-04, 0.2767314E-04, 0.3277707E-04, &
0.3864321E-04, 0.4547839E-04, 0.5328839E-04, 0.6234301E-04, &
0.7263268E-04, 0.8450696E-04, 0.9793231E-04, 0.1133587E-03, &
0.1307170E-03, 0.1505832E-03, 0.1728373E-03, 0.1982122E-03, &
0.2266389E-03, 0.2592220E-03, 0.2957792E-03, 0.3376068E-03, &
0.3844381E-03, 0.4379281E-03, 0.4976965E-03, 0.5658476E-03, &
0.6418494E-03, 0.7287094E-03, 0.8261995E-03, 0.9380076E-03, &
0.1063498E-02, 0.1207423E-02, 0.1369594E-02, 0.1557141E-02, &
0.1769657E-02, 0.2015887E-02, 0.2295520E-02, 0.2620143E-02, &
0.2989651E-02, 0.3419469E-02, 0.3909867E-02, 0.4481491E-02, &
0.5135272E-02, 0.5898971E-02, 0.6774619E-02, 0.7799763E-02/
DATA P2/ &
0.8978218E-02, 0.1036103E-01, 0.1195488E-01, 0.1382957E-01, &
0.1599631E-01, 0.1855114E-01, 0.2151235E-01, 0.2501293E-01, &
0.2908220E-01, 0.3390544E-01, 0.3952926E-01, 0.4621349E-01, &
0.5403168E-01, 0.6330472E-01, 0.7406807E-01, 0.8677983E-01, &
0.1015345E+00, 0.1189603E+00, 0.1391863E+00, 0.1630739E+00, &
0.1908004E+00, 0.2235461E+00, 0.2609410E+00, 0.3036404E+00, &
0.3513750E+00, 0.4055375E+00, 0.4656677E+00, 0.5335132E+00, &
0.6083618E+00, 0.6923932E+00, 0.7845676E+00, 0.8875882E+00, &
0.1000000E+01/
DATA O3HI1/ &
.55,.50,.45,.45,.40,.35,.35,.30,.30,.30, &
.55,.51,.46,.47,.42,.38,.37,.36,.35,.35, &
.55,.53,.48,.49,.44,.42,.41,.40,.38,.38, &
.60,.55,.52,.52,.50,.47,.46,.44,.42,.41, &
.65,.60,.55,.56,.53,.52,.50,.48,.45,.45, &
.75,.65,.60,.60,.55,.55,.55,.50,.48,.47, &
.80,.75,.75,.75,.70,.70,.65,.63,.60,.60, &
.90,.85,.85,.80,.80,.75,.75,.74,.72,.71, &
1.10,1.05,1.00,.90,.90,.90,.85,.83,.80,.80, &
1.40,1.30,1.25,1.25,1.25,1.20,1.15,1.10,1.05,1.00, &
1.7,1.7,1.6,1.6,1.6,1.6,1.6,1.6,1.5,1.5, &
2.1,2.0,1.9,1.9,1.9,1.8,1.8,1.8,1.7,1.7, &
2.4,2.3,2.2,2.2,2.2,2.1,2.1,2.1,2.0,2.0, &
2.7,2.5,2.5,2.5,2.5,2.5,2.4,2.4,2.3,2.3, &
2.9,2.8,2.7