source: codes/icosagcm/trunk/src/etat0_dcmip4.f90 @ 342

Last change on this file since 342 was 286, checked in by dubos, 10 years ago

Partial etat0 cleanup (removed calls to xyz2lonlat)
File size: 7.6 KB
Line 
1MODULE etat0_dcmip4_mod
2  USE icosa
3  PRIVATE
4  REAL(rstd),PARAMETER :: eta0=0.252
5  REAL(rstd),PARAMETER :: etat=0.2
6  REAL(rstd),PARAMETER :: ps0=1e5
7  REAL(rstd),PARAMETER :: u0=35
8  REAL(rstd),PARAMETER :: T0=288
9  REAL(rstd),PARAMETER :: DeltaT=4.8e5
10  REAL(rstd),PARAMETER :: Rd=287
11  REAL(rstd),PARAMETER :: Gamma=0.005
12  REAL(rstd),PARAMETER :: up0=1
13  REAL(rstd) :: latw=2*Pi/9
14!$OMP THREADPRIVATE(latw)
15  REAL(rstd) :: pw=34000
16!$OMP THREADPRIVATE(pw)
17  REAL(rstd) :: q0=0.021
18!$OMP THREADPRIVATE(q0)
19  REAL(rstd) :: lonc
20!$OMP THREADPRIVATE(lonc)
21  REAL(rstd) :: latc
22!$OMP THREADPRIVATE(latc)
23
24  INTEGER,SAVE :: testcase
25!$OMP THREADPRIVATE(testcase) 
26
27  PUBLIC  etat0
28CONTAINS
29 
30   
31   
32  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
33  USE icosa
34  IMPLICIT NONE
35    TYPE(t_field),POINTER :: f_ps(:)
36    TYPE(t_field),POINTER :: f_phis(:)
37    TYPE(t_field),POINTER :: f_theta_rhodz(:)
38    TYPE(t_field),POINTER :: f_u(:)
39    TYPE(t_field),POINTER :: f_q(:)
40 
41    REAL(rstd),POINTER :: ps(:)
42    REAL(rstd),POINTER :: phis(:)
43    REAL(rstd),POINTER :: theta_rhodz(:,:)
44    REAL(rstd),POINTER :: u(:,:)
45    REAL(rstd),POINTER :: q(:,:,:)
46    INTEGER :: ind
47
48    testcase=1
49    CALL getin("dcmip4_testcase",testcase)
50   
51    DO ind=1,ndomain
52      IF (.NOT. assigned_domain(ind)) CYCLE
53      CALL swap_dimensions(ind)
54      CALL swap_geometry(ind)
55      ps=f_ps(ind)
56      phis=f_phis(ind)
57      theta_rhodz=f_theta_rhodz(ind)
58      u=f_u(ind)
59      q=f_q(ind)
60      CALL compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q)
61    ENDDO
62
63  END SUBROUTINE etat0
64 
65  SUBROUTINE compute_etat0_dcmip4(ps, phis, theta_rhodz, u, q)
66  USE icosa
67  USE disvert_mod
68  USE pression_mod
69  USE exner_mod
70  USE geopotential_mod
71  USE theta2theta_rhodz_mod
72  IMPLICIT NONE 
73  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
74  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
75  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
76  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
77  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
78 
79  INTEGER :: i,j,l,ij
80  REAL(rstd) :: theta(iim*jjm,llm)
81  REAL(rstd) :: Y(iim*jjm,llm)
82  REAL(rstd) :: vort
83  REAL(rstd) :: eta(llm)
84  REAL(rstd) :: etav(llm)
85  REAL(rstd) :: etas, etavs
86  REAL(rstd) :: lon,lat
87  REAL(rstd) :: ulon(3)
88  REAL(rstd) :: ep(3), norm_ep
89  REAL(rstd) :: Tave, T
90  REAL(rstd) :: phis_ave
91  REAL(rstd) :: V0(3)
92  REAL(rstd) :: r2
93  REAL(rstd) :: utot
94  REAL(rstd) :: lonx,latx
95  REAL(rstd) :: dthetaodeta_ave, dthetaodeta, dthetaodlat, duodeta, K, r
96 
97    lonc=Pi/9
98    latc=2*Pi/9 
99 
100    DO l=1,llm
101      eta(l)= 0.5 *( ap(l)/preff+bp(l) + ap(l+1)/preff+bp(l+1) )
102      etav(l)=(eta(l)-eta0)*Pi/2
103    ENDDO
104    etas=ap(1)/preff+bp(1)
105    etavs=(etas-eta0)*Pi/2
106
107    DO j=jj_begin,jj_end
108      DO i=ii_begin,ii_end
109        ij=(j-1)*iim+i
110        ps(ij)=ps0
111      ENDDO
112    ENDDO
113   
114   
115    CALL lonlat2xyz(lonc,latc,V0)
116
117    u(:,:)=1e10     
118    DO l=1,llm
119      DO j=jj_begin-1,jj_end+1
120        DO i=ii_begin-1,ii_end+1
121          ij=(j-1)*iim+i
122         
123          lon=lon_e(ij+u_right) ; lat=lat_e(ij+u_right)
124          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
125          r=radius*acos(K)
126          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2)
127          u(ij+u_right,l) = utot * sum(elon_e(ij+u_right,:) * ep_e(ij+u_right,:))
128
129          lon=lon_e(ij+u_lup) ; lat=lat_e(ij+u_lup)
130          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
131          r=radius*acos(K)
132          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2)
133          u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:))
134
135          lon=lon_e(ij+u_ldown) ; lat=lat_e(ij+u_ldown)
136          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
137          r=radius*acos(K)
138          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2)
139          u(ij+u_ldown,l) = utot * sum(elon_e(ij+u_ldown,:) * ep_e(ij+u_ldown,:))
140
141       ENDDO
142     ENDDO
143    ENDDO 
144     
145     
146     DO l=1,llm
147       Tave=T0*eta(l)**(Rd*Gamma/g)
148       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
149       DO j=jj_begin-1,jj_end+1
150         DO i=ii_begin-1,ii_end+1
151           ij=(j-1)*iim+i
152           lat=lat_i(ij)
153            Y(ij,l)=((-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5     &
154                                + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega)
155            T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l)
156           
157            theta(ij,l)=T*eta(l)**(-kappa)
158
159          ENDDO
160       ENDDO
161     ENDDO
162     
163     
164     phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
165     DO j=jj_begin,jj_end
166       DO i=ii_begin,ii_end
167         ij=(j-1)*iim+i
168         lat=lat_i(ij)
169         phis(ij)=phis_ave+u0*cos(etavs)**1.5*( (-2*sin(lat)**6 * (cos(lat)**2+1./3) + 10./63 )*u0*cos(etavs)**1.5  &
170                                           +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega )
171!         phis(ij)=phis_ave+u0*cos(etavs)**1.5
172
173       ENDDO
174     ENDDO
175
176    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
177
178
179    IF (testcase==1) THEN
180   
181      q(:,:,1)=theta(:,:)
182      IF(nqtot>2) q(:,:,3)=1.
183
184      DO l=1,llm
185         dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1)
186         IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa)  &
187                                                          + kappa * (etat-eta(l))**5 * eta(l)**(-kappa-1))
188         DO j=jj_begin,jj_end
189           DO i=ii_begin,ii_end
190             ij=(j-1)*iim+i
191             lon=lon_i(ij) ; lat=lat_i(ij)
192             dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) & 
193                                        + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l)                    & 
194                                   - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)&
195                                   - 9./8.  * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))   &
196                                                                                  * (-2*sin(lat)**6*(cos(lat)**2+1./3.)+10./63.) 
197             dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5                                   &
198                         *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) &
199                        + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat)))
200                         
201             duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l))
202           
203             K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
204             r=radius*acos(K)
205             vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2)                  & 
206                    + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) &
207                                                            -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2)))
208             q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta))
209             IF(nqtot>3) q(ij,l,4)=cos(lon)*cos(lat)
210           ENDDO
211         ENDDO
212       ENDDO       
213
214    ELSE IF (testcase==2) THEN
215   
216      DO l=1,llm
217         DO j=jj_begin,jj_end
218           DO i=ii_begin,ii_end
219             ij=(j-1)*iim+i
220             q(ij,l,1)=q0*exp(-(lat_i(ij)/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2)
221           ENDDO
222         ENDDO
223       ENDDO
224   
225    ENDIF       
226
227       
228  END SUBROUTINE compute_etat0_dcmip4
229 
230END MODULE etat0_dcmip4_mod
Note: See TracBrowser for help on using the repository browser.