source: codes/icosagcm/trunk/src/initial/etat0_dcmip4.f90

Last change on this file was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

Fixed iso c bindings
fixed warnings with -Wall -Wno-aliasing -Wno-unused -Wno-unused-dummy-argument -Wno-maybe-uninitialized -Wno-tabs warnings
Removed all unused variables (-Wunused-variable)
vector%dot_product is now dot_product_3d to avoid compilation warning "dot_product shadows intrinsic" with GCC
File size: 4.7 KB
Line 
1MODULE etat0_dcmip4_mod
2  USE icosa
3  IMPLICIT NONE
4  PRIVATE
5 
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
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 
19  INTEGER,SAVE :: testcase
20  !$OMP THREADPRIVATE(testcase) 
21 
22  PUBLIC getin_etat0, compute_etat0
23
24CONTAINS
25
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
34    testcase=1
35    CALL getin("dcmip4_testcase",testcase)
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
41    USE omp_para
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)
51   
52    INTEGER :: l,ij
53    REAL(rstd) :: etal, etavl, etas, etavs, sinlat, coslat, &
54         Y, Tave, T, phis_ave, vort, utot, &
55         dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r
56   
57    etas=ap(1)/preff+bp(1)
58    etavs=(etas-eta0)*Pi/2
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
67    ENDDO
68   
69    DO l=ll_begin,ll_end
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))
84         
85          K=sin(latc)*sinlat+cos(latc)*coslat*cos(lon(ij)-lonc)
86          r=radius*acos(K)
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
119   
120  END SUBROUTINE compute_etat0
121 
122END MODULE etat0_dcmip4_mod
Note: See TracBrowser for help on using the repository browser.