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

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

Partial etat0 cleanup (removed calls to xyz2lonlat)
File size: 8.5 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 dist_lonlat(lon0,lat0,lon_i(n),lat_i(n),rr1)  ! GC distance from center
191            rr1 = radius*rr1
192            IF ( rr1 .LT. R0 ) then
193              hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
194            ELSE
195              hx(n,l)=0.0
196            END IF
197          END DO
198        END DO
199      END DO
200   
201    END SUBROUTINE cosine_bell_1
202
203!==============================================================================
204
205    SUBROUTINE   cosine_bell_2(hx) 
206    IMPLICIT NONE       
207    REAL(rstd) :: hx(iim*jjm,llm)
208   
209      DO l=1,llm
210        DO j=jj_begin-1,jj_end+1
211          DO i=ii_begin-1,ii_end+1
212            n=(j-1)*iim+i
213            CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center
214            rr1 = radius*rr1
215            CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center
216            rr2 = radius*rr2
217            dd1t1 = rr1/rt
218            dd1t1 = dd1t1*dd1t1
219            dd1t2 = (zrl(l) - zc)/zt
220            dd1t2 = dd1t2*dd1t2 
221            dd1 = dd1t1 + dd1t2 
222            dd1 = Min(1.0,dd1) 
223            dd2t1 = rr2/rt
224            dd2t1 = dd2t1*dd2t1
225            dd2 = dd2t1 + dd1t2 
226            dd2 = Min(1.0,dd2) 
227            hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
228          END DO
229        END DO
230      END DO
231
232    END SUBROUTINE cosine_bell_2
233
234!=============================================================================
235
236    SUBROUTINE slotted_cylinders(hx) 
237    IMPLICIT NONE
238    REAL(rstd) :: hx(iim*jjm,llm)
239   
240      DO l=1,llm
241        DO j=jj_begin-1,jj_end+1
242          DO i=ii_begin-1,ii_end+1
243            n=(j-1)*iim+i
244            CALL dist_lonlat(lonc1,latc1,lon_i(n),lat_i(n),rr1)  ! GC distance from center
245            rr1 = radius*rr1
246            CALL dist_lonlat(lonc2,latc2,lon_i(n),lat_i(n),rr2)  ! GC distance from center
247            rr2 = radius*rr2
248            dd1t1 = rr1/rt
249            dd1t1 = dd1t1*dd1t1
250            dd1t2 = (zrl(l) - zc)/zt
251            dd1t2 = dd1t2*dd1t2 
252            dd1 = dd1t1 + dd1t2 
253            dd2t1 = rr2/rt
254            dd2t1 = dd2t1*dd2t1
255            dd2 = dd2t1 + dd1t2
256           
257            IF ( dd1 .LT. 0.5 ) Then
258              hx(n,l) = 1.0
259            ELSEIF ( dd2 .LT. 0.5 ) Then
260              hx(n,l) = 1.0
261            ELSE
262              hx(n,l) = 0.1 
263            END IF 
264   
265            IF ( zrl(l) .GT. zc ) Then
266              IF ( ABS(latc1 - lat) .LT. 0.125 ) Then
267                hx(n,l)= 0.1
268              ENDIF
269            ENDIF 
270 
271          ENDDO
272       END DO
273      END DO
274
275    END SUBROUTINE slotted_cylinders 
276
277!==============================================================================
278
279    SUBROUTINE hadleyq(hx) 
280      IMPLICIT NONE
281      REAL(rstd)::hx(iim*jjm,llm) 
282      REAL(rstd),PARAMETER:: zz1=2000.,zz2=5000.,zz0=0.5*(zz1+zz2)
283     
284      DO l=1,llm
285         IF ( ( zz1 .LT. zrl(l) ) .and. ( zrl(l) .LT. zz2 ) )  THEN
286            hx(:,l) = 0.5*(1. + cos(2*pi*(zrl(l)-zz0)/(zz2-zz1))) 
287         ELSE
288            hx(:,l) = 0.0 
289         END IF
290      END DO
291    END SUBROUTINE hadleyq
292
293  END SUBROUTINE compute_etat0_ncar
294
295END MODULE etat0_dcmip1_mod
Note: See TracBrowser for help on using the repository browser.