[78] | 1 | MODULE etat0_dcmip4_mod |
---|
[62] | 2 | USE icosa |
---|
[346] | 3 | IMPLICIT NONE |
---|
[62] | 4 | PRIVATE |
---|
[346] | 5 | |
---|
[62] | 6 | REAL(rstd),PARAMETER :: eta0=0.252 |
---|
| 7 | REAL(rstd),PARAMETER :: etat=0.2 |
---|
| 8 | REAL(rstd),PARAMETER :: ps0=1e5 |
---|
| 9 | REAL(rstd),PARAMETER :: u0=35 |
---|
| 10 | REAL(rstd),PARAMETER :: T0=288 |
---|
| 11 | REAL(rstd),PARAMETER :: DeltaT=4.8e5 |
---|
| 12 | REAL(rstd),PARAMETER :: Rd=287 |
---|
| 13 | REAL(rstd),PARAMETER :: Gamma=0.005 |
---|
| 14 | REAL(rstd),PARAMETER :: up0=1 |
---|
[346] | 15 | REAL(rstd),PARAMETER :: lonc=Pi/9, latc=2*Pi/9, latw=2*Pi/9 |
---|
| 16 | REAL(rstd),PARAMETER :: pw=34000 |
---|
| 17 | REAL(rstd),PARAMETER :: q0=0.021 |
---|
| 18 | |
---|
[186] | 19 | INTEGER,SAVE :: testcase |
---|
[346] | 20 | !$OMP THREADPRIVATE(testcase) |
---|
| 21 | |
---|
| 22 | PUBLIC getin_etat0, compute_etat0 |
---|
[186] | 23 | |
---|
[62] | 24 | CONTAINS |
---|
[186] | 25 | |
---|
[346] | 26 | SUBROUTINE getin_etat0 |
---|
| 27 | USE mpipara, ONLY : is_mpi_root |
---|
| 28 | IF(nqtot<2) THEN |
---|
| 29 | IF (is_mpi_root) THEN |
---|
| 30 | PRINT *, "nqtot must be at least 2 for test case DCMIP4" |
---|
| 31 | END IF |
---|
| 32 | STOP |
---|
| 33 | END IF |
---|
[186] | 34 | testcase=1 |
---|
| 35 | CALL getin("dcmip4_testcase",testcase) |
---|
[346] | 36 | END SUBROUTINE getin_etat0 |
---|
| 37 | |
---|
| 38 | SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q) |
---|
| 39 | USE icosa |
---|
| 40 | USE disvert_mod |
---|
[353] | 41 | USE omp_para |
---|
[346] | 42 | INTEGER, INTENT(IN) :: ngrid |
---|
| 43 | REAL(rstd),INTENT(IN) :: lon(ngrid) |
---|
| 44 | REAL(rstd),INTENT(IN) :: lat(ngrid) |
---|
| 45 | REAL(rstd),INTENT(OUT) :: phis(ngrid) |
---|
| 46 | REAL(rstd),INTENT(OUT) :: ps(ngrid) |
---|
| 47 | REAL(rstd),INTENT(OUT) :: temp(ngrid,llm) |
---|
| 48 | REAL(rstd),INTENT(OUT) :: ulon(ngrid,llm) |
---|
| 49 | REAL(rstd),INTENT(OUT) :: ulat(ngrid,llm) |
---|
| 50 | REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot) |
---|
[62] | 51 | |
---|
[346] | 52 | INTEGER :: l,ij |
---|
| 53 | REAL(rstd) :: etal, etavl, etas, etavs, sinlat, coslat, & |
---|
[899] | 54 | Y, Tave, T, phis_ave, vort, utot, & |
---|
[346] | 55 | dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r |
---|
| 56 | |
---|
[67] | 57 | etas=ap(1)/preff+bp(1) |
---|
[62] | 58 | etavs=(etas-eta0)*Pi/2 |
---|
[346] | 59 | phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g)) |
---|
| 60 | |
---|
| 61 | DO ij=1,ngrid |
---|
| 62 | sinlat=SIN(lat(ij)) |
---|
| 63 | coslat=COS(lat(ij)) |
---|
| 64 | phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sinlat**6 * (coslat**2+1./3) + 10./63 )*u0*cos(etavs)**1.5 & |
---|
| 65 | +(8./5*coslat**3 * (sinlat**2 + 2./3) - Pi/4)*radius*Omega ) |
---|
| 66 | ps(ij)=ps0 |
---|
[62] | 67 | ENDDO |
---|
| 68 | |
---|
[353] | 69 | DO l=ll_begin,ll_end |
---|
[346] | 70 | etal = 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) ) |
---|
| 71 | etavl=(etal-eta0)*Pi/2 |
---|
| 72 | |
---|
| 73 | Tave=T0*etal**(Rd*Gamma/g) |
---|
| 74 | dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* etal**(Rd*Gamma/g-kappa-1) |
---|
| 75 | IF (etat>etal) THEN |
---|
| 76 | Tave=Tave+DeltaT*(etat-etal)**5 |
---|
| 77 | dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-etal)**4 * etal**(-kappa) & |
---|
| 78 | + kappa * (etat-etal)**5 * etal**(-kappa-1)) |
---|
| 79 | END IF |
---|
| 80 | |
---|
| 81 | DO ij=1,ngrid |
---|
| 82 | sinlat=SIN(lat(ij)) |
---|
| 83 | coslat=COS(lat(ij)) |
---|
[62] | 84 | |
---|
[346] | 85 | K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc) |
---|
[62] | 86 | r=radius*acos(K) |
---|
[346] | 87 | utot=u0*cos(etavl)**1.5*sin(2*lat(ij))**2 + up0*exp(-(r/(0.1*radius))**2) |
---|
| 88 | ulon(ij,l) = utot |
---|
| 89 | ulat(ij,l) = 0. |
---|
| 90 | Y = ((-2*sinlat**6*(coslat**2+1./3)+10./63)*2*u0*cos(etavl)**1.5 & |
---|
| 91 | + (8./5*coslat**3*(sinlat**2+2./3)-Pi/4)*radius*Omega) |
---|
| 92 | T = Tave + 0.75*(etal*Pi*u0/Rd)*sin(etavl)*cos(etavl)**0.5 * Y |
---|
| 93 | temp(ij,l)=T |
---|
| 94 | |
---|
| 95 | IF (testcase==1) THEN |
---|
| 96 | q(ij,l,1)=T*etal**(-kappa) |
---|
| 97 | IF(nqtot>2) q(ij,l,3)=1. |
---|
| 98 | dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*etal**(-kappa)*sin(etavl)*cos(etavl)**0.5 * Y & |
---|
| 99 | + 3/8. * Pi**2*u0/Rd * etal**(1-kappa) * cos(etavl)**1.5 * Y & |
---|
| 100 | - 3./16. * Pi**2 * u0 /Rd * etal**(1-kappa) * sin(etavl)**2 * cos(etavl)**(-0.5) *Y & |
---|
| 101 | - 9./8. * Pi**2 * u0 /Rd * etal**(1-kappa) * sin(etavl)**2 * cos(etavl) & |
---|
| 102 | * (-2*sinlat**6*(coslat**2+1./3.)+10./63.) |
---|
| 103 | dthetaodlat=3./4.*Pi*u0/Rd*etal**(1-kappa)*sin(etavl)*cos(etavl)**0.5 & |
---|
| 104 | *( 2*u0*cos(etavl)**1.5 * ( -12 * coslat*sinlat**5*(coslat**2+1./3.)+4*coslat*sinlat**7) & |
---|
| 105 | + radius*omega*(-24./5. * sinlat * coslat**2 * (sinlat**2 + 2./3.) + 16./5. * coslat**4 * sinlat)) |
---|
| 106 | |
---|
| 107 | duodeta=-u0 * sin(2*lat(ij))**2 * 3./4.*Pi * cos(etavl)**0.5 * sin(etavl) |
---|
| 108 | |
---|
| 109 | vort = -4*u0/radius*cos(etavl)**1.5 * sinlat * coslat * (2.-5.*sinlat**2) & |
---|
| 110 | + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat(ij))-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*coslat & |
---|
| 111 | -cos(latc)*sinlat*cos(lon(ij)-lonc))/(sqrt(1-K**2))) |
---|
| 112 | q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sinlat*omega+vort)*dthetaodeta)) |
---|
| 113 | IF(nqtot>3) q(ij,l,4)=cos(lon(ij))*coslat |
---|
| 114 | ELSE IF (testcase==2) THEN |
---|
| 115 | q(ij,l,1)=q0*exp(-(lat(ij)/latw)**4)*exp(-((etal-1)*preff/pw)**2) |
---|
| 116 | END IF |
---|
| 117 | END DO |
---|
| 118 | END DO |
---|
[62] | 119 | |
---|
[346] | 120 | END SUBROUTINE compute_etat0 |
---|
[62] | 121 | |
---|
[78] | 122 | END MODULE etat0_dcmip4_mod |
---|