MODULE m_Parameters IMPLICIT NONE PRIVATE PUBLIC :: Get_Parameters INTEGER, PARAMETER, PUBLIC :: ik4=6, rk4=6, ik8=14, rk8=14 INTEGER, PARAMETER, PUBLIC :: ki4=SELECTED_INT_KIND(ik4), kr4=SELECTED_REAL_KIND(rk4) INTEGER, PARAMETER, PUBLIC :: ki8=SELECTED_INT_KIND(ik8), kr8=SELECTED_REAL_KIND(rk8) INTEGER (KIND=ki8), PUBLIC :: k, Kmax, Kmaxm, Kmaxp REAL (KIND=kr8), PUBLIC :: rd, gv, cp, ps, pt, ts, tt, zs, zt, rk, dpt, gvbyrd, tiso REAL (KIND=kr8) :: pa(1), ta(1), za(1) REAL (KIND=kr8), DIMENSION (:), ALLOCATABLE, PUBLIC :: PrsLev, prsint, delprs, t0, z0, gama CHARACTER (LEN=5), PUBLIC :: csvm CHARACTER (LEN=25), PUBLIC :: dnprs CHARACTER (LEN=25), PUBLIC :: dnfnb CHARACTER (LEN=29), PUBLIC :: dnfng CHARACTER (LEN=29), PUBLIC :: dnfns CHARACTER (LEN=128), PUBLIC :: dirprs, dirvmd LOGICAL (KIND=ki4) :: isotermica NAMELIST /vmdnml/ csvm, Kmax, isotermica, tiso, dirprs, dirvmd CONTAINS SUBROUTINE Get_Parameters IMPLICIT NONE ! Default Values for Namelist csvm='ERA_5' Kmax=37 isotermica=.FALSE. tiso=287.15_kr8 dirprs='../../Pressure_Levels/ ' dirvmd='../dataout/ ' ! Other Parameters ts=288.15_kr8 zs=0.0_kr8 tiso=300.0_kr8 rd=287.05_kr8 cp=1005.0_kr8 gv=9.80665_kr8 rk=rd/cp gvbyrd=1000.0_kr8*gv/rd OPEN (UNIT=55, FILE='./NameList', STATUS='OLD') READ (UNIT=55, NML=vmdnml) CLOSE (UNIT=55) Kmaxm = Kmax-1 Kmaxp = Kmax+1 dnfnb='Vertical_Modes_'//csvm//'.L ' dnfng='Vertical_Functions_'//csvm//'.L ' dnfns='Vertical_Standards_'//csvm//'.L ' dnprs='PressureLevels_'//csvm//'.L ' WRITE (UNIT=dnfnb(23:25), FMT='(I3.3)') Kmax WRITE (UNIT=dnfng(27:29), FMT='(I3.3)') Kmax WRITE (UNIT=dnfns(27:29), FMT='(I3.3)') Kmax WRITE (UNIT=dnprs(23:25), FMT='(I3.3)') Kmax WRITE (UNIT=*, FMT='(/,A)') ' $vmdnml' WRITE (UNIT=*, FMT='(A,I5)') ' Kmax = ', Kmax WRITE (UNIT=*, FMT='(A,L5)') ' isotermica = ', isotermica WRITE (UNIT=*, FMT='(A,F12.5)') ' tiso = ',tiso WRITE (UNIT=*, FMT='(A)') ' dirprs = '//TRIM(dirprs) WRITE (UNIT=*, FMT='(A)') ' dirvmd = '//TRIM(dirvmd) WRITE (UNIT=*, FMT='(A)') ' /' WRITE (UNIT=*, FMT='(A)') ' dnfnb = '//dnfnb WRITE (UNIT=*, FMT='(A)') ' dnfng = '//dnfng WRITE (UNIT=*, FMT='(A)') ' dnfns = '//dnfns WRITE (UNIT=*, FMT='(A)') ' dnprs = '//dnprs WRITE (UNIT=*, FMT='(/,A,I5)') ' ik4 = ', ik4 WRITE (UNIT=*, FMT='(A,I5)') ' rk4 = ', rk4 WRITE (UNIT=*, FMT='(A,I5)') ' ki4 = ', ki4 WRITE (UNIT=*, FMT='(A,I5)') ' kr4 = ', kr4 WRITE (UNIT=*, FMT='(/,A,I5)') ' ik8 = ', ik8 WRITE (UNIT=*, FMT='(A,I5)') ' rk8 = ', rk8 WRITE (UNIT=*, FMT='(A,I5)') ' ki8 = ', ki8 WRITE (UNIT=*, FMT='(A,I5)') ' kr8 = ', kr8 WRITE (UNIT=*, FMT='(/,A)') ' ' WRITE (UNIT=*, FMT='(A,F9.3)') ' rd = ', rd WRITE (UNIT=*, FMT='(A,F9.3)') ' cp = ', cp WRITE (UNIT=*, FMT='(A,F9.3)') ' gv = ', gv WRITE (UNIT=*, FMT='(A,F9.3)') ' rk = ', rk WRITE (UNIT=*, FMT='(A,F9.3)') ' gr = ', gvbyrd WRITE (UNIT=*, FMT='(/,A)') ' ' ALLOCATE (PrsLev(Kmaxp), prsint(Kmaxp), delprs(Kmaxp)) ALLOCATE (t0(Kmaxp), z0(Kmaxp), gama(Kmaxp)) ! Reading Pressure Levels ! Pressure: First Surface, Last Top ! Revert Levels in Reading: Top to Surface OPEN (UNIT=90, FILE=TRIM(dirprs)//dnprs, STATUS='OLD') READ (UNIT=90, FMT='(10F9.3)') PrsLev(Kmax:1:-1) READ (UNIT=90, FMT='(F9.3)') ps CLOSE (UNIT=90) ! Top Pressure and Pressure of the Lower Extra Level Below Surface dpt=0.5_kr8*(PrsLev(2)-PrsLev(1)) IF (PrsLev(1) > dpt) THEN pt=PrsLev(1)-dpt ELSE pt=0.5_kr8*PrsLev(1) ENDIF PrsLev(Kmaxp)=ps+(ps-PrsLev(Kmax)) ! Delta Pressure on Layers delprs(1)=2.0_kr8*(PrsLev(1)-pt) DO k=2,Kmax delprs(k)=PrsLev(k)-PrsLev(k-1) ENDDO delprs(Kmaxp)=2.0_kr8*(ps-PrsLev(Kmax)) ! Pressure on Interfaces prsint(1)=pt DO k=2,Kmax prsint(k)=0.5_kr8*(PrsLev(k-1)+PrsLev(k)) ENDDO prsint(Kmaxp)=ps IF (isotermica) THEN zs=0.0_kr8 DO k=1,Kmaxp z0(k)=zs+((rd*tiso)/gv)*LOG(ps/PrsLev(k)) END DO t0=tiso zt=zs+((rd*tiso)/gv)*LOG(ps/pt) tt=tiso ELSE ! Surface Temperature and Geopotential Height pa(1)=ps CALL Standard_Atmosphere (pa, ta, za, 1_ki8) ts=ta(1) zs=za(1) ! Top Temperature and Geopotential Height pa(1)=pt CALL Standard_Atmosphere (pa, ta, za, 1_ki8) tt=ta(1) zt=za(1) ! Compute Basic State Temperature and Geopotential Height on Layers CALL Standard_Atmosphere (PrsLev, t0, z0, Kmaxp) ENDIF ! Compute Static Stability on Interfaces gama(1)=(rd/(pt*gv))*((rd/(pt*cp))*tt-(t0(1)-tt)/(0.5_kr8*delprs(1))) DO k=2,Kmax gama(k)=(rd/(prsint(k)*gv))*((rd/(prsint(k)*cp))*0.5_kr8*(t0(k-1)+t0(k))-(t0(k)-t0(k-1))/delprs(k)) enddo gama(Kmaxp)=(rd/(ps*gv))*((rd/(ps*cp))*ts-(ts-t0(Kmax))/(0.5_kr8*delprs(Kmaxp))) WRITE (UNIT=*, FMT='(/,A,F9.3)') ' ps = ',ps WRITE (UNIT=*, FMT='(A,F9.3)') ' zs = ',zs WRITE (UNIT=*, FMT='(A,F9.3,/)') ' ts = ',ts WRITE (UNIT=*, FMT='(A,F9.3)') ' pt = ',pt WRITE (UNIT=*, FMT='(A,F9.3)') ' zt = ',zt WRITE (UNIT=*, FMT='(A,F9.3)') ' tt = ',tt WRITE (UNIT=*, FMT='(/,A)') ' Pressure on Layers:' WRITE (UNIT=*, FMT='(8F11.3)') PrsLev WRITE (UNIT=*, FMT='(A)') ' ' WRITE (UNIT=*, FMT='(/,A)') ' Pressure on Interfaces:' WRITE (UNIT=*, FMT='(8F11.3)') prsint WRITE (UNIT=*, FMT='(A)') ' ' WRITE (UNIT=*, FMT='(/,A)') ' Delta Pressure on Layers:' WRITE (UNIT=*, FMT='(8F11.3)') delprs WRITE (UNIT=*, FMT='(/,A)') ' Temperature on Layers:' WRITE (UNIT=*, FMT='(8F11.3)') t0 WRITE (UNIT=*, FMT='(/,A)') ' Geopotencial Height on Layers:' WRITE (UNIT=*, FMT='(8F11.3)') z0 WRITE (UNIT=*, FMT='(/,A)') ' Static Stability on Interfaces:' WRITE (UNIT=*, FMT='(1P5E12.5)') gama WRITE (UNIT=*, FMT='(A)') ' ' !IF (PrsLev(1) /= 0.0_kr8) STOP END SUBROUTINE Get_Parameters SUBROUTINE Standard_Atmosphere (p, t, z, np) IMPLICIT NONE ! Computes the Temperature And Geopotential Height of the 1976 ! U.S. Standard Atmosphere, up to 86 km of Geometric Altitude, ! Given the Pressure Level. INTEGER (KIND=ki8), INTENT(IN) :: np REAL (KIND=kr8), INTENT(IN) :: p(np) REAL (KIND=kr8), INTENT(OUT) :: t(np) REAL (KIND=kr8), INTENT(OUT) :: z(np) ! Arguments: ! np (Input): Number of Requested Pressure Levels ! p (Input): Pressure of Requested LevelS (hPa) ! t (Output): Temperature of Requested Levels (K) ! z (Output): Geopotential Height of Requested Levels (m) ! Local: ! nstand: Number of Significant Levels of the Standard Atmosphere INTEGER (KIND=ki8), PARAMETER :: nstand=8 ! Counters INTEGER (KIND=ki8) :: n, i, j, k ! Auxiliary Constante REAL (KIND=kr8) :: cte ! Surface Temperature (K) and Pressure (hPa) and the ! Aceleration of Gravity Divided by Dry Air Gas Constant (K/km) !JPB REAL (KIND=kr8) :: ts=288.15, ps=1000.0, gvbyrd=34.1631947 ! Standard Geopotential Heights (km) (84.852 Geop km = 86 Geometric km) ! Convertion of Geometric Altitude (alt) Into Geopotential Height (h): ! h = alt*rearth/(alt+rearth) , rearth = Earth Mean Radius REAL (KIND=kr8), DIMENSION (nstand) :: & hstand=(/0.0, 11.0, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852/) ! Standard Vertical Virtual Temperature Gradients (K/km) ! Note that gstand is Minus the Lapse-Rate (= - dt/dh) REAL (KIND=kr8), DIMENSION (nstand-1) :: gstand=(/-6.5, 0.0, 1.0, 2.8, 0.0, -2.8, -2.0/) ! Temperature and Pressure of the Significant Levels REAL (KIND=kr8) :: tstand(nstand), pstand(nstand) tstand(1)=ts DO n=2,nstand tstand(n)=tstand(n-1)+gstand(n-1)*(hstand(n)-hstand(n-1)) END DO pstand(1)=ps DO n=2,nstand IF (n == 3 .OR. n == 6) THEN pstand(n)=pstand(n-1)*EXP(-(gvbyrd/tstand(n-1))* & (hstand(n)-hstand(n-1))) ELSE pstand(n)=pstand(n-1)*EXP(-(gvbyrd/gstand(n-1))* & LOG(tstand(n)/tstand(n-1))) END IF END DO DO n=1,np i=1 j=nstand 10 k=(i+j)/2 IF (p(n) > pstand(k)) THEN j=k ELSE i=k END IF IF (j <= i+1) GO TO 20 GO TO 10 20 cte=gstand(i)/gvbyrd IF (cte == 0.0) THEN t(n)=tstand(i) z(n)=1000.0*(hstand(i)-(tstand(i)/gvbyrd)*LOG(p(n)/pstand(i))) ELSE t(n)=tstand(i)*EXP(-cte*LOG(p(n)/pstand(i))) z(n)=1000.0*(hstand(i)+(t(n)-tstand(i))/gstand(i)) END IF END DO END SUBROUTINE Standard_Atmosphere END MODULE m_Parameters