source: codes/icosagcm/trunk/src/etat0_dcmip5.f90 @ 176

Last change on this file since 176 was 114, checked in by ymipsl, 12 years ago

correction for initial temperature on dcmip testcase 5

YM
File size: 9.2 KB
Line 
1MODULE etat0_dcmip5_mod
2  USE icosa
3  PRIVATE
4  REAL(rstd),PARAMETER :: zt=15000
5  REAL(rstd),PARAMETER :: q0=0.021
6  REAL(rstd),PARAMETER :: qt=1e-11
7  REAL(rstd),PARAMETER :: T0=302.15
8  REAL(rstd),PARAMETER :: Ts=302.15
9  REAL(rstd),PARAMETER :: zq1=3000
10  REAL(rstd),PARAMETER :: zq2=8000
11  REAL(rstd),PARAMETER :: Gamma=0.007
12  REAL(rstd),PARAMETER :: pb=101500
13  REAL(rstd),PARAMETER :: Deltap=1115
14  REAL(rstd),PARAMETER :: rp=282000
15  REAL(rstd),PARAMETER :: zp=7000
16  REAL(rstd),PARAMETER :: epsilon=1e-25
17  REAL(rstd),PARAMETER :: Rd=287
18  REAL(rstd) :: lonc
19  REAL(rstd) :: latc
20  TYPE(t_field),POINTER :: f_out_i(:)
21  REAL(rstd),POINTER :: out_i(:,:)
22
23  PUBLIC  etat0
24CONTAINS
25 
26   
27   
28  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
29  USE icosa
30  IMPLICIT NONE
31    TYPE(t_field),POINTER :: f_ps(:)
32    TYPE(t_field),POINTER :: f_phis(:)
33    TYPE(t_field),POINTER :: f_theta_rhodz(:)
34    TYPE(t_field),POINTER :: f_u(:)
35    TYPE(t_field),POINTER :: f_q(:)
36 
37    REAL(rstd),POINTER :: ps(:)
38    REAL(rstd),POINTER :: phis(:)
39    REAL(rstd),POINTER :: theta_rhodz(:,:)
40    REAL(rstd),POINTER :: u(:,:)
41    REAL(rstd),POINTER :: q(:,:,:)
42    INTEGER :: ind
43 
44    CALL allocate_field(f_out_i,field_t,type_real,llm)
45 
46    DO ind=1,ndomain
47      CALL swap_dimensions(ind)
48      CALL swap_geometry(ind)
49      ps=f_ps(ind)
50      phis=f_phis(ind)
51      theta_rhodz=f_theta_rhodz(ind)
52      u=f_u(ind)
53      q=f_q(ind)
54      out_i=f_out_i(ind)
55      CALL compute_etat0_dcmip5(ps, phis, theta_rhodz, u, q)
56    ENDDO
57    CALL writefield("out_i",f_out_i)
58    CALL deallocate_field(f_out_i)
59  END SUBROUTINE etat0
60 
61  SUBROUTINE compute_etat0_dcmip5(ps, phis, theta_rhodz, ue, q)
62  USE icosa
63  USE pression_mod
64  USE wind_mod
65  USE theta2theta_rhodz_mod
66  USE exner_mod
67  IMPLICIT NONE 
68  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
69  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
70  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
71  REAL(rstd),INTENT(OUT) :: ue(3*iim*jjm,llm)
72  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
73 
74  INTEGER :: i,j,l,ij
75  INTEGER :: testcase
76  REAL(rstd) :: Tv0, Tvt 
77  REAL(rstd) :: r 
78  REAL(rstd) :: z(iim*jjm,llm),zz 
79  REAL(rstd) :: p(iim*jjm,llm+1) 
80  REAL(rstd) :: Tave,Tp 
81  REAL(rstd) :: T(iim*jjm,llm) 
82  REAL(rstd) :: vt 
83  REAL(rstd) :: d1,d2,d 
84  REAL(rstd) :: ulon(3*iim*jjm,llm) 
85  REAL(rstd) :: ulat(3*iim*jjm,llm)
86  REAL(rstd) :: lon,lat
87  REAL(rstd) :: pks(iim*jjm) 
88  REAL(rstd) :: pk(iim*jjm,llm) 
89
90
91    latc=Pi/18
92    lonc=Pi 
93
94    Tv0=T0/(1+0.608*q0)
95    Tvt=Tv0-Gamma*zt
96 
97    testcase=1
98    CALL getin("dcmip5_testcase",testcase)
99   
100    phis(:)=0.
101   
102    DO j=jj_begin,jj_end
103      DO i=ii_begin,ii_end
104        ij=(j-1)*iim+i
105        CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
106        r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
107        ps(ij)=pb-DeltaP*exp(-(r/rp)**1.5)
108      ENDDO
109    ENDDO
110   
111    CALL compute_pression(ps,p,0)
112
113    DO l=1,llm
114      DO j=jj_begin,jj_end
115        DO i=ii_begin,ii_end
116          ij=(j-1)*iim+i
117          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
118          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
119          CALL compute_z(0.5*(p(ij,l)+p(ij,l+1)),ps(ij),r,z(ij,l))
120        ENDDO
121      ENDDO
122    ENDDO
123    out_i=z
124   
125    DO l=1,llm
126      DO j=jj_begin-1,jj_end+1
127        DO i=ii_begin-1,ii_end+1
128          ij=(j-1)*iim+i
129          IF (z(ij,l)<=zt) THEN
130            q(ij,l,1)=q0*exp(-z(ij,l)/zq1)*exp(-(z(ij,l)/zq2)**2)
131          ELSE
132            q(ij,l,1)=qt
133          ENDIF
134        ENDDO
135      ENDDO
136    ENDDO
137
138   
139     DO l=1,llm
140      DO j=jj_begin-1,jj_end+1
141        DO i=ii_begin-1,ii_end+1
142          ij=(j-1)*iim+i
143          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
144          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
145          IF (z(ij,l)<=zt) THEN
146            Tave=Tv0-Gamma*z(ij,l)
147            Tp=(Tv0-Gamma*z(ij,l)) * ( ( 1 + 2*Rd*(Tv0-Gamma*z(ij,l))*z(ij,l) /                           &
148                                            (g*zp**2*(1-pb/Deltap*exp((r/rp)**1.5+(z(ij,l)/zp)**2))))**(-1) - 1) 
149          ELSE
150            Tave=Tvt
151            Tp=0
152          ENDIF
153          T(ij,l)=Tave+Tp 
154          out_i(ij,l)=Tave+Tp
155        ENDDO
156      ENDDO
157    ENDDO
158
159    DO l=1,llm
160      DO j=jj_begin-1,jj_end+1
161        DO i=ii_begin-1,ii_end+1
162          ij=(j-1)*iim+i
163
164          CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
165          zz=0.5*(z(ij,l)+z(ij+t_right,l))
166          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
167          IF (zz<=zt) THEN
168            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
169                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
170                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
171          ELSE
172            vt = 0
173          ENDIF
174          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
175          d2=cos(latc)*sin(lon-lonc)
176          d=MAX(1e-25,sqrt(d1**2+d2**2))
177          ulon(ij+u_right,l)=vt*d1/d
178          ulat(ij+u_right,l)=vt*d2/d
179
180
181         
182          CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
183          zz=0.5*(z(ij,l)+z(ij+t_lup,l))
184          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
185          IF (zz<=zt) THEN
186            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
187                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
188                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
189
190         ELSE
191            vt = 0
192          ENDIF
193          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
194          d2=cos(latc)*sin(lon-lonc)
195          d=MAX(1e-25,sqrt(d1**2+d2**2))
196          ulon(ij+u_lup,l)=vt*d1/d
197          ulat(ij+u_lup,l)=vt*d2/d
198         
199         
200         
201         
202          CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
203          zz=0.5*(z(ij,l)+z(ij+t_ldown,l))
204          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
205          IF (zz<=zt) THEN
206            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
207                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
208                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
209
210          ELSE
211            vt = 0
212          ENDIF
213          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
214          d2=cos(latc)*sin(lon-lonc)
215          d=MAX(1e-25,sqrt(d1**2+d2**2))
216          ulon(ij+u_ldown,l)=vt*d1/d
217          ulat(ij+u_ldown,l)=vt*d2/d
218         
219
220
221          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
222          zz=z(ij,l)
223          r=radius*acos(sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc))
224          IF (zz<=zt) THEN
225            vt=-omega*sin(latc)*r+sqrt( (omega*sin(latc)*r)**2 -  1.5*(r/rp)**1.5 * (Tv0-Gamma*zz)*Rd    &
226                                                                 / (1 + 2*Rd*(Tv0-Gamma*zz)*zz/(g*zp**2)        &
227                                                                      - pb/Deltap * exp((r/rp)**1.5 + (zz/zp)**2 )))
228          ELSE
229            vt = 0
230          ENDIF
231          d1=sin(latc)*cos(lat)-cos(latc)*sin(lat)*cos(lon-lonc)
232          d2=cos(latc)*sin(lon-lonc)
233          d=MAX(1e-25,sqrt(d1**2+d2**2))
234!          ulon(ij+u_ldown,l)=vt*d1/d
235!          ulat(ij+u_ldown,l)=vt*d2/d
236!          out_i(ij,l)=sqrt((vt*d1/d)**2+(vt*d2/d)**2)
237        ENDDO
238      ENDDO
239    ENDDO
240   
241    CALL compute_wind_perp_from_lonlat_compound(ulon, ulat, ue)
242    CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
243
244  END SUBROUTINE compute_etat0_dcmip5
245 
246  SUBROUTINE compute_z(pmodel,ps,r,z)
247  USE icosa
248  IMPLICIT NONE
249     REAL(rstd),PARAMETER   :: epsilon0=2e-13
250     REAL(rstd),INTENT(IN)  :: pmodel
251     REAL(rstd),INTENT(IN)  :: ps
252     REAL(rstd),INTENT(IN)  :: r
253     REAL(rstd),INTENT(OUT) :: z
254     
255     REAL(rstd) :: Tv0,Tvt,pt
256     REAL(rstd) :: pave, pp, dpdz
257     REAL(rstd) :: znp1
258     REAL(rstd) :: epsilon
259     REAL(rstd) :: p
260     INTEGER  :: n
261         
262     Tv0=T0/(1+0.608*q0)
263     Tvt=Tv0-Gamma*zt
264     pt=pb*(Tvt/Tv0)**(g/(Rd*Gamma))
265
266     IF (pmodel>pt) THEN
267       z=Tv0/Gamma*(1-(pmodel/ps)**(Rd*Gamma/g))
268     ELSE
269       z=zt+Rd*Tvt/g*log(Pt/pmodel)
270     ENDIF
271     
272!     PRINT*,pmodel,pt,z,zt
273     IF (z>zt .OR. r>1000000) return
274     
275     epsilon=1
276     n=0
277     
278     DO WHILE(epsilon>epsilon0 .AND. n<20)
279 
280       IF (z<zt) THEN
281         pave=pb*((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))
282         pp= -Deltap * exp(-(r/rp)**1.5-(z/zp)**2) *  ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma))
283       ELSE
284         pave=pt*exp(g*(zt-z)/(Rd*Tvt))
285         pp=0
286       ENDIF
287       p=pave+pp
288     
289       dpdz =  2*Deltap*z/zp**2 * exp(-(r/rp)**1.5) * exp(-(z/zp)**2) * ((Tv0-Gamma*z)/Tv0)**(g/(Rd*Gamma))  &
290             - g/(Rd*Tv0) * (preff - Deltap*exp(-(r/rp)**1.5) * exp(-(z/zp)**2)) * ( (Tv0-Gamma*z)/Tv0 )**(g/(Rd*Gamma)-1)
291           
292       znp1 = z - ( pmodel - p ) / (-dpdz)
293       
294       epsilon=ABS( (znp1-z)/znp1 )
295!       PRINT *,"epsilon",epsilon,r,z,znp1
296!       PRINT *,"----",n,pmodel,ps,pave,dpdz
297       z=znp1
298       n=n+1
299     ENDDO 
300   
301   END SUBROUTINE compute_z
302
303END MODULE etat0_dcmip5_mod
Note: See TracBrowser for help on using the repository browser.