source: codes/icosagcm/trunk/src/etat0_dcmip1.f90 @ 195

Last change on this file since 195 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: 8.6 KB
Line 
1MODULE etat0_dcmip1_mod
2  USE icosa
3  PRIVATE
4
5  REAL(rstd), SAVE  :: h0=1.
6!$OMP THREADPRIVATE(h0)
7  REAL(rstd), SAVE  :: lon0=3*pi/2
8!$OMP THREADPRIVATE(lon0)
9  REAL(rstd), SAVE  :: lat0=0.0
10!$OMP THREADPRIVATE(lat0)
11  REAL(rstd), SAVE  :: alpha=0.0
12!$OMP THREADPRIVATE(alpha)
13  REAL(rstd), SAVE  :: R0 
14!$OMP THREADPRIVATE(R0)
15  REAL(rstd), SAVE  :: lat1=0.
16!$OMP THREADPRIVATE(lat1)
17  REAL(rstd), SAVE  :: lat2=0.
18!$OMP THREADPRIVATE(lat2)
19  REAL(rstd), SAVE  :: lon1=pi/6
20!$OMP THREADPRIVATE(lon1)
21  REAL(rstd), SAVE  :: lon2=-pi/6
22!$OMP THREADPRIVATE(lon2)
23  REAL(rstd), SAVE  :: latc1=0.
24!$OMP THREADPRIVATE(latc1)
25  REAL(rstd), SAVE  :: latc2=0.
26!$OMP THREADPRIVATE(latc2)
27  REAL(rstd), SAVE  :: lonc1=5*pi/6
28!$OMP THREADPRIVATE(lonc1)
29  REAL(rstd), SAVE  :: lonc2=7*pi/6
30!$OMP THREADPRIVATE(lonc2)
31  REAL(rstd), SAVE  :: zt=1000.0
32!$OMP THREADPRIVATE(zt)
33  REAL(rstd), SAVE  :: rt
34!$OMP THREADPRIVATE(rt)
35  REAL(rstd), SAVE  :: zc=5000.0
36!$OMP THREADPRIVATE(zc)
37
38  PUBLIC etat0
39 
40CONTAINS
41 
42  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
43  USE icosa
44  IMPLICIT NONE
45    TYPE(t_field),POINTER :: f_ps(:)
46    TYPE(t_field),POINTER :: f_phis(:)
47    TYPE(t_field),POINTER :: f_theta_rhodz(:)
48    TYPE(t_field),POINTER :: f_u(:)
49    TYPE(t_field),POINTER :: f_q(:)
50 
51    REAL(rstd),POINTER :: ps(:)
52    REAL(rstd),POINTER :: phis(:)
53    REAL(rstd),POINTER :: theta_rhodz(:,:)
54    REAL(rstd),POINTER :: u(:,:)
55    REAL(rstd),POINTER :: q(:,:,:)
56    CHARACTER(len=255) :: dcmip1_adv_shape
57    INTEGER :: ind
58   
59    R0=radius*0.5
60    rt=radius*0.5
61    dcmip1_adv_shape='cos_bell'
62    CALL getin('dcmip1_shape',dcmip1_adv_shape)
63   
64    DO ind=1,ndomain
65      IF (.NOT. assigned_domain(ind)) CYCLE
66      CALL swap_dimensions(ind)
67      CALL swap_geometry(ind)
68      ps=f_ps(ind)
69      phis=f_phis(ind)
70      theta_rhodz=f_theta_rhodz(ind)
71      u=f_u(ind)
72      q=f_q(ind)
73
74
75      SELECT CASE(TRIM(dcmip1_adv_shape))
76      CASE('const')
77         CALL compute_etat0_ncar(1,ps, phis, theta_rhodz, u, q(:,:,1))
78      CASE('cos_bell')
79         CALL compute_etat0_ncar(2,ps, phis, theta_rhodz, u, q(:,:,1))
80      CASE('slotted_cyl')
81         CALL compute_etat0_ncar(3,ps, phis, theta_rhodz, u, q(:,:,1))
82      CASE('dbl_cos_bell_q1')
83         CALL compute_etat0_ncar(4,ps, phis, theta_rhodz, u, q(:,:,1))
84      CASE('dbl_cos_bell_q2')
85         CALL compute_etat0_ncar(5,ps, phis, theta_rhodz, u, q(:,:,1))
86      CASE('complement')
87         CALL compute_etat0_ncar(6,ps, phis, theta_rhodz, u, q(:,:,1))
88      CASE('hadley')  ! hadley like meridional circulation
89         CALL compute_etat0_ncar(7,ps, phis, theta_rhodz, u, q(:,:,1))
90      CASE('dcmip11')
91         IF(nqtot==5) THEN
92            CALL compute_etat0_ncar(4,ps, phis, theta_rhodz, u, q(:,:,1))
93            CALL compute_etat0_ncar(5,ps, phis, theta_rhodz, u, q(:,:,2))
94            CALL compute_etat0_ncar(3,ps, phis, theta_rhodz, u, q(:,:,3))
95            CALL compute_etat0_ncar(6,ps, phis, theta_rhodz, u, q(:,:,4))
96            CALL compute_etat0_ncar(1,ps, phis, theta_rhodz, u, q(:,:,5))
97         ELSE
98            PRINT *,'Error : etat0_dcmip=dcmip11 and nqtot = ',nqtot,' .'
99            PRINT *,'nqtot must be equal to 5 when etat0_dcmip=dcmip11'
100            STOP
101         END IF
102      CASE DEFAULT
103         PRINT *, 'Bad selector for variable dcmip1_adv_shape : <', TRIM(dcmip1_adv_shape),   &
104              '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
105              '<dbl_cos_bell_q2>, <complement>, <hadley>'
106         STOP
107      END SELECT
108
109    ENDDO
110
111  END SUBROUTINE etat0
112 
113  SUBROUTINE compute_etat0_ncar(icase, ps, phis, theta_rhodz, u, q)
114  USE icosa
115  USE disvert_mod
116  USE pression_mod
117  USE exner_mod
118  USE geopotential_mod
119  USE theta2theta_rhodz_mod
120  IMPLICIT NONE 
121  INTEGER, INTENT(in) :: icase
122  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
123  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
124  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
125  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
126  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm) 
127
128  REAL(rstd) :: qxt1(iim*jjm,llm)
129  REAL(rstd) :: lon, lat
130  REAL(rstd) ::dd1,dd2,dd1t1,dd1t2,dd2t1
131  REAL(rstd) :: pr, zr(llm+1), zrl(llm)
132  REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
133  REAL(rstd) :: X2(3),X1(3)
134  INTEGER :: i,j,n,l
135
136  u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0
137 
138  DO l=1, llm+1
139     pr = ap(l) + bp(l)*ncar_p0
140     zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0) 
141  ENDDO
142 
143  DO l=1, llm 
144     zrl(l) = 0.5*(zr(l) + zr(l+1)) 
145  END DO
146 
147  SELECT CASE(icase)
148  CASE(1)
149     q=1
150  CASE(2)
151     !--------------------------------------------- SINGLE COSINE BELL
152     CALL cosine_bell_1(q)     
153  CASE(3)
154     CALL slotted_cylinders(q) 
155  CASE(4)
156     PRINT *, 'Double cosine bell'
157     CALL cosine_bell_2(q)     
158  CASE(5)
159     CALL cosine_bell_2(q) 
160     DO l=1,llm
161        q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l)
162     END DO     
163  CASE(6)
164     ! tracer such that, in combination with the other tracer fields
165     ! with weight (3/10), the sum is equal to one
166     CALL cosine_bell_2(qxt1)
167     DO l = 1,llm
168        q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l) 
169     END DO
170     q = q + qxt1 
171     CALL slotted_cylinders(qxt1)
172     q = q + qxt1 
173     q = 1. - q*0.3
174  CASE(7)  ! hadley like meridional circulation
175     CALL hadleyq(q) 
176  END SELECT
177     
178CONTAINS
179
180!======================================================================
181
182    SUBROUTINE cosine_bell_1(hx)
183    IMPLICIT NONE
184    REAL(rstd) :: hx(iim*jjm,llm)
185   
186      DO l=1,llm 
187        DO j=jj_begin-1,jj_end+1
188          DO i=ii_begin-1,ii_end+1
189            n=(j-1)*iim+i
190            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
191            CALL dist_lonlat(lon0,lat0,lon,lat,rr1)  ! GC distance from center
192            rr1 = radius*rr1
193            IF ( rr1 .LT. R0 ) then
194              hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
195            ELSE
196              hx(n,l)=0.0
197            END IF
198          END DO
199        END DO
200      END DO
201   
202    END SUBROUTINE cosine_bell_1
203
204!==============================================================================
205
206    SUBROUTINE   cosine_bell_2(hx) 
207    IMPLICIT NONE       
208    REAL(rstd) :: hx(iim*jjm,llm)
209   
210      DO l=1,llm
211        DO j=jj_begin-1,jj_end+1
212          DO i=ii_begin-1,ii_end+1
213            n=(j-1)*iim+i
214            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
215            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center
216            rr1 = radius*rr1
217            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center
218            rr2 = radius*rr2
219            dd1t1 = rr1/rt
220            dd1t1 = dd1t1*dd1t1
221            dd1t2 = (zrl(l) - zc)/zt
222            dd1t2 = dd1t2*dd1t2 
223            dd1 = dd1t1 + dd1t2 
224            dd1 = Min(1.0,dd1) 
225            dd2t1 = rr2/rt
226            dd2t1 = dd2t1*dd2t1
227            dd2 = dd2t1 + dd1t2 
228            dd2 = Min(1.0,dd2) 
229            hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
230          END DO
231        END DO
232      END DO
233
234    END SUBROUTINE cosine_bell_2
235
236!=============================================================================
237
238    SUBROUTINE slotted_cylinders(hx) 
239    IMPLICIT NONE
240    REAL(rstd) :: hx(iim*jjm,llm)
241   
242      DO l=1,llm
243        DO j=jj_begin-1,jj_end+1
244          DO i=ii_begin-1,ii_end+1
245            n=(j-1)*iim+i
246            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
247            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center
248            rr1 = radius*rr1
249            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center
250            rr2 = radius*rr2
251            dd1t1 = rr1/rt
252            dd1t1 = dd1t1*dd1t1
253            dd1t2 = (zrl(l) - zc)/zt
254            dd1t2 = dd1t2*dd1t2 
255            dd1 = dd1t1 + dd1t2 
256            dd2t1 = rr2/rt
257            dd2t1 = dd2t1*dd2t1
258            dd2 = dd2t1 + dd1t2
259           
260            IF ( dd1 .LT. 0.5 ) Then
261              hx(n,l) = 1.0
262            ELSEIF ( dd2 .LT. 0.5 ) Then
263              hx(n,l) = 1.0
264            ELSE
265              hx(n,l) = 0.1 
266            END IF 
267   
268            IF ( zrl(l) .GT. zc ) Then
269              IF ( ABS(latc1 - lat) .LT. 0.125 ) Then
270                hx(n,l)= 0.1
271              ENDIF
272            ENDIF 
273 
274          ENDDO
275       END DO
276      END DO
277
278    END SUBROUTINE slotted_cylinders 
279
280!==============================================================================
281
282    SUBROUTINE hadleyq(hx) 
283    IMPLICIT NONE
284    REAL(rstd)::hx(iim*jjm,llm) 
285    REAL(rstd),PARAMETER:: zz1=3500.,zz2=6500.,zz0=0.5*(zz1+zz2)
286   
287      DO l=1,llm
288        IF ( ( zz1 .LT. zrl(l) ) .and. ( zrl(l) .LT. zz2 ) )  THEN
289          hx(:,l) = 0.5*(1. + cos(0.002*pi*(zrl(l)-zz0)/3.)) 
290        ELSE
291          hx(:,l) = 0.0 
292        END IF
293      END DO
294    END SUBROUTINE hadleyq 
295
296  END SUBROUTINE compute_etat0_ncar
297 
298
299END MODULE etat0_dcmip1_mod
Note: See TracBrowser for help on using the repository browser.