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

Last change on this file since 186 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM
File size: 7.8 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          CALL xyz2lonlat(xyz_e(ij+u_right,:)/radius,lon,lat)
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
130          CALL xyz2lonlat(xyz_e(ij+u_lup,:)/radius,lon,lat)
131          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
132          r=radius*acos(K)
133          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2)
134          u(ij+u_lup,l) = utot * sum(elon_e(ij+u_lup,:) * ep_e(ij+u_lup,:))
135
136          CALL xyz2lonlat(xyz_e(ij+u_ldown,:)/radius,lon,lat)
137          K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
138          r=radius*acos(K)
139          utot=u0*cos(etav(l))**1.5*sin(2*lat)**2 + up0*exp(-(r/(0.1*radius))**2)
140          u(ij+u_ldown,l) = utot * sum(elon_e(ij+u_ldown,:) * ep_e(ij+u_ldown,:))
141
142       ENDDO
143     ENDDO
144    ENDDO 
145     
146     
147     DO l=1,llm
148       Tave=T0*eta(l)**(Rd*Gamma/g)
149       IF (etat>eta(l)) Tave=Tave+DeltaT*(etat-eta(l))**5
150       DO j=jj_begin-1,jj_end+1
151         DO i=ii_begin-1,ii_end+1
152           ij=(j-1)*iim+i
153           CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
154           
155            Y(ij,l)=((-2*sin(lat)**6*(cos(lat)**2+1./3)+10./63)*2*u0*cos(etav(l))**1.5     &
156                                + (8./5*cos(lat)**3*(sin(lat)**2+2./3)-Pi/4)*radius*Omega)
157            T=Tave+ 0.75*(eta(l)*Pi*u0/Rd)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l)
158           
159            theta(ij,l)=T*eta(l)**(-kappa)
160
161          ENDDO
162       ENDDO
163     ENDDO
164     
165     
166     phis_ave=T0*g/Gamma*(1-etas**(Rd*Gamma/g))
167     DO j=jj_begin,jj_end
168       DO i=ii_begin,ii_end
169         ij=(j-1)*iim+i
170         CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
171         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  &
172                                           +(8./5*cos(lat)**3 * (sin(lat)**2 + 2./3) - Pi/4)*radius*Omega )
173!         phis(ij)=phis_ave+u0*cos(etavs)**1.5
174
175       ENDDO
176     ENDDO
177
178    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,0)
179
180
181    IF (testcase==1) THEN
182   
183      q(:,:,1)=theta(:,:)
184      IF(nqtot>2) q(:,:,3)=1.
185
186      DO l=1,llm
187         dthetaodeta_ave = T0 *( Rd*Gamma/g - kappa)* eta(l)**(Rd*Gamma/g-kappa-1)
188         IF (etat>eta(l)) dthetaodeta_ave = dthetaodeta_ave - DeltaT * ( 5*(etat-eta(l))**4 * eta(l)**(-kappa)  &
189                                                          + kappa * (etat-eta(l))**5 * eta(l)**(-kappa-1))
190         DO j=jj_begin,jj_end
191           DO i=ii_begin,ii_end
192             ij=(j-1)*iim+i
193             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
194             dthetaodeta=dthetaodeta_ave + 3./4. * Pi * u0/Rd*(1-kappa)*eta(l)**(-kappa)*sin(etav(l))*cos(etav(l))**0.5 * Y(ij,l) & 
195                                        + 3/8. * Pi**2*u0/Rd * eta(l)**(1-kappa) * cos(etav(l))**1.5 * Y(ij,l)                    & 
196                                   - 3./16. * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))**(-0.5) *Y(ij,l)&
197                                   - 9./8.  * Pi**2 * u0 /Rd * eta(l)**(1-kappa) * sin(etav(l))**2 * cos(etav(l))   &
198                                                                                  * (-2*sin(lat)**6*(cos(lat)**2+1./3.)+10./63.) 
199             dthetaodlat=3./4.*Pi*u0/Rd*eta(l)**(1-kappa)*sin(etav(l))*cos(etav(l))**0.5                                   &
200                         *( 2*u0*cos(etav(l))**1.5 * ( -12 * cos(lat)*sin(lat)**5*(cos(lat)**2+1./3.)+4*cos(lat)*sin(lat)**7) &
201                        + radius*omega*(-24./5. * sin(lat) * cos(lat)**2 * (sin(lat)**2 + 2./3.) + 16./5. * cos(lat)**4 * sin(lat)))
202                         
203             duodeta=-u0 * sin(2*lat)**2 * 3./4.*Pi * cos(etav(l))**0.5 * sin(etav(l))
204           
205             K=sin(latc)*sin(lat)+cos(latc)*cos(lat)*cos(lon-lonc)
206             r=radius*acos(K)
207             vort = -4*u0/radius*cos(etav(l))**1.5 * sin(lat) * cos(lat) * (2.-5.*sin(lat)**2)                  & 
208                    + up0/radius*exp(-(r/(0.1*radius))**2) * (tan(lat)-2*(radius/(0.1*radius))**2 * acos(K) * (sin(latc)*cos(lat) &
209                                                            -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2)))
210             q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta))
211             IF(nqtot>3) q(ij,l,4)=cos(lon)*cos(lat)
212           ENDDO
213         ENDDO
214       ENDDO       
215
216    ELSE IF (testcase==2) THEN
217   
218      DO l=1,llm
219         DO j=jj_begin,jj_end
220           DO i=ii_begin,ii_end
221             ij=(j-1)*iim+i
222             CALL xyz2lonlat(xyz_i(ij,:)/radius,lon,lat)
223             q(ij,l,1)=q0*exp(-(lat/latw)**4)*exp(-((eta(l)-1)*preff/pw)**2)
224           ENDDO
225         ENDDO
226       ENDDO
227   
228    ENDIF       
229
230       
231  END SUBROUTINE compute_etat0_dcmip4
232 
233END MODULE etat0_dcmip4_mod
Note: See TracBrowser for help on using the repository browser.