source: codes/icosagcm/trunk/src/etat0_ncar.f90 @ 20

Last change on this file since 20 was 19, checked in by ymipsl, 12 years ago

Simplify the management of the module.

YM

File size: 6.7 KB
Line 
1MODULE etat0_ncar_mod
2  USE icosa
3  PRIVATE
4
5  REAL(rstd), PARAMETER :: h0=1.
6  REAL(rstd), PARAMETER :: lon0=3*pi/2
7  REAL(rstd), PARAMETER :: lat0=0.0
8  REAL(rstd), PARAMETER :: alpha=0.0
9  REAL(rstd), PARAMETER :: R0 = radius*0.5
10  REAL(rstd), PARAMETER :: lat1=0.
11  REAL(rstd), PARAMETER :: lat2=0.
12  REAL(rstd), PARAMETER :: lon1=pi/6
13  REAL(rstd), PARAMETER :: lon2=-pi/6
14  REAL(rstd), PARAMETER :: latc1=0.
15  REAL(rstd), PARAMETER :: latc2=0.
16  REAL(rstd), PARAMETER :: lonc1=5*pi/6
17  REAL(rstd), PARAMETER :: lonc2=7*pi/6
18  REAL(rstd), PARAMETER :: zt=1000.0
19  REAL(rstd), PARAMETER :: rt=radius*0.5
20  REAL(rstd), PARAMETER :: zc=5000.0
21  REAL(rstd), PARAMETER :: ps0=100000.0
22  REAL(rstd), PARAMETER :: T0=300
23  REAL(rstd), PARAMETER :: R_d=287.0   
24
25  PUBLIC etat0
26 
27CONTAINS
28 
29  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
30  USE icosa
31  IMPLICIT NONE
32    TYPE(t_field),POINTER :: f_ps(:)
33    TYPE(t_field),POINTER :: f_phis(:)
34    TYPE(t_field),POINTER :: f_theta_rhodz(:)
35    TYPE(t_field),POINTER :: f_u(:)
36    TYPE(t_field),POINTER :: f_q(:)
37 
38    REAL(rstd),POINTER :: ps(:)
39    REAL(rstd),POINTER :: phis(:)
40    REAL(rstd),POINTER :: theta_rhodz(:,:)
41    REAL(rstd),POINTER :: u(:,:)
42    REAL(rstd),POINTER :: q(:,:,:)
43    INTEGER :: ind
44   
45    DO ind=1,ndomain
46      CALL swap_dimensions(ind)
47      CALL swap_geometry(ind)
48      ps=f_ps(ind)
49      phis=f_phis(ind)
50      theta_rhodz=f_theta_rhodz(ind)
51      u=f_u(ind)
52      q=f_q(ind)
53      CALL compute_etat0_ncar(ps, phis, theta_rhodz, u, q(:,:,1))
54    ENDDO
55
56  END SUBROUTINE etat0
57 
58  SUBROUTINE compute_etat0_ncar(ps, phis, theta_rhodz, u, q)
59  USE icosa
60  USE disvert_mod
61  USE pression_mod
62  USE exner_mod
63  USE geopotential_mod
64  USE theta2theta_rhodz_mod
65  IMPLICIT NONE 
66  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
67  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
68  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
69  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
70  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm) 
71
72
73  REAL(rstd) :: qxt1(iim*jjm,llm)
74  REAL(rstd) :: lon, lat
75  REAL(rstd) ::dd1,dd2,dd1t1,dd1t2,dd2t1
76  REAL(rstd) :: pr, zr(llm+1), zrl(llm)
77  REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
78  REAL(rstd) :: X2(3),X1(3)
79  INTEGER :: i,j,n,l
80  INTEGER :: testcase
81
82    u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ps0
83       
84    DO l=1, llm+1
85      pr = ap(l) + bp(l)*ps0
86      zr(l) = -R_d*T0/g*log(pr/ps0) 
87    ENDDO   
88
89    DO l=1, llm 
90      zrl(l) = 0.5*(zr(l) + zr(l+1)) 
91    END DO
92
93    testcase=5
94    CALL getin('ncar_test_case',testcase)
95 
96    SELECT CASE(testcase)
97!--------------------------------------------- SINGLE COSINE BELL
98      CASE(0)
99        CALL cosine_bell_1(q) 
100       
101      CASE(1)
102        CALL cosine_bell_2(q) 
103       
104      CASE(2)
105        CALL cosine_bell_2(q) 
106        DO l=1,llm
107          q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l)
108        END DO 
109       
110      CASE(3)
111        CALL slotted_cylinders(q)       
112
113      CASE(4) 
114        CALL cosine_bell_2(qxt1)
115        DO l = 1,llm
116          q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l) 
117        END DO
118        q = q + qxt1 
119        CALL slotted_cylinders(qxt1)
120        q = q + qxt1 
121        q = 1. - q*0.3 
122         
123      CASE(5)  ! hadley like meridional circulation
124        CALL hadleyq(q) 
125       
126      CASE DEFAULT
127        PRINT*,"no such initial profile"
128        STOP
129       
130      END SELECT
131 
132  CONTAINS
133
134!======================================================================
135
136    SUBROUTINE cosine_bell_1(hx)
137    IMPLICIT NONE
138    REAL(rstd) :: hx(iim*jjm,llm)
139   
140      DO l=1,llm 
141        DO j=jj_begin-1,jj_end+1
142          DO i=ii_begin-1,ii_end+1
143            n=(j-1)*iim+i
144            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
145            CALL dist_lonlat(lon0,lat0,lon,lat,rr1)  ! GC distance from center
146            rr1 = radius*rr1
147            IF ( rr1 .LT. R0 ) then
148              hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
149            ELSE
150              hx(n,l)=0.0
151            END IF
152          END DO
153        END DO
154      END DO
155   
156    END SUBROUTINE cosine_bell_1
157
158!==============================================================================
159
160    SUBROUTINE   cosine_bell_2(hx) 
161    IMPLICIT NONE       
162    REAL(rstd) :: hx(iim*jjm,llm)
163   
164      DO l=1,llm
165        DO j=jj_begin-1,jj_end+1
166          DO i=ii_begin-1,ii_end+1
167            n=(j-1)*iim+i
168            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
169            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center
170            rr1 = radius*rr1
171            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center
172            rr2 = radius*rr2
173            dd1t1 = rr1/rt
174            dd1t1 = dd1t1*dd1t1
175            dd1t2 = (zrl(l) - zc)/zt
176            dd1t2 = dd1t2*dd1t2 
177            dd1 = dd1t1 + dd1t2 
178            dd1 = Min(1.0,dd1) 
179            dd2t1 = rr2/rt
180            dd2t1 = dd2t1*dd2t1
181            dd2 = dd2t1 + dd1t2 
182            dd2 = Min(1.0,dd2) 
183            hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
184          END DO
185        END DO
186      END DO
187
188    END SUBROUTINE cosine_bell_2
189
190!=============================================================================
191
192    SUBROUTINE slotted_cylinders(hx) 
193    IMPLICIT NONE
194    REAL(rstd) :: hx(iim*jjm,llm)
195   
196      DO l=1,llm
197        DO j=jj_begin-1,jj_end+1
198          DO i=ii_begin-1,ii_end+1
199            n=(j-1)*iim+i
200            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
201            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center
202            rr1 = radius*rr1
203            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center
204            rr2 = radius*rr2
205            dd1t1 = rr1/rt
206            dd1t1 = dd1t1*dd1t1
207            dd1t2 = (zrl(l) - zc)/zt
208            dd1t2 = dd1t2*dd1t2 
209            dd1 = dd1t1 + dd1t2 
210            dd2t1 = rr2/rt
211            dd2t1 = dd2t1*dd2t1
212            dd2 = dd2t1 + dd1t2
213           
214            IF ( dd1 .LT. 0.5 ) Then
215              hx(n,l) = 1.0
216            ELSEIF ( dd2 .LT. 0.5 ) Then
217              hx(n,l) = 1.0
218            ELSE
219              hx(n,l) = 0.1 
220            END IF 
221   
222            IF ( zrl(l) .GT. zc ) Then
223              IF ( ABS(latc1 - lat) .LT. 0.125 ) Then
224                hx(n,l)= 0.0
225              ENDIF
226            ENDIF 
227 
228          ENDDO
229       END DO
230      END DO
231
232    END SUBROUTINE slotted_cylinders 
233
234!==============================================================================
235
236    SUBROUTINE hadleyq(hx) 
237    IMPLICIT NONE
238    REAL(rstd)::hx(iim*jjm,llm) 
239    REAL(rstd),PARAMETER:: zz1=3500.,zz2=6500.,zz0=0.5*(zz1+zz2)
240   
241      DO l=1,llm
242        IF ( ( zz1 .LT. zrl(l) ) .and. ( zrl(l) .LT. zz2 ) )  THEN
243          hx(:,l) = 0.5*(1. + cos(0.002*pi*(zrl(l)-zz0)/3.)) 
244        ELSE
245          hx(:,l) = 0.0 
246        END IF
247      END DO
248    END SUBROUTINE hadleyq 
249
250  END SUBROUTINE compute_etat0_ncar
251 
252
253END MODULE etat0_ncar_mod
Note: See TracBrowser for help on using the repository browser.