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

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

Some correction due to externalisation of earth const parameters

YM

File size: 6.9 KB
Line 
1MODULE etat0_ncar_mod
2  USE icosa
3  PRIVATE
4
5  REAL(rstd), SAVE  :: h0=1.
6  REAL(rstd), SAVE  :: lon0=3*pi/2
7  REAL(rstd), SAVE  :: lat0=0.0
8  REAL(rstd), SAVE  :: alpha=0.0
9  REAL(rstd), SAVE  :: R0 
10  REAL(rstd), SAVE  :: lat1=0.
11  REAL(rstd), SAVE  :: lat2=0.
12  REAL(rstd), SAVE  :: lon1=pi/6
13  REAL(rstd), SAVE  :: lon2=-pi/6
14  REAL(rstd), SAVE  :: latc1=0.
15  REAL(rstd), SAVE  :: latc2=0.
16  REAL(rstd), SAVE  :: lonc1=5*pi/6
17  REAL(rstd), SAVE  :: lonc2=7*pi/6
18  REAL(rstd), SAVE  :: zt=1000.0
19  REAL(rstd), SAVE  :: rt
20  REAL(rstd), SAVE  :: zc=5000.0
21
22  PUBLIC etat0
23 
24CONTAINS
25 
26  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
27  USE icosa
28  IMPLICIT NONE
29    TYPE(t_field),POINTER :: f_ps(:)
30    TYPE(t_field),POINTER :: f_phis(:)
31    TYPE(t_field),POINTER :: f_theta_rhodz(:)
32    TYPE(t_field),POINTER :: f_u(:)
33    TYPE(t_field),POINTER :: f_q(:)
34 
35    REAL(rstd),POINTER :: ps(:)
36    REAL(rstd),POINTER :: phis(:)
37    REAL(rstd),POINTER :: theta_rhodz(:,:)
38    REAL(rstd),POINTER :: u(:,:)
39    REAL(rstd),POINTER :: q(:,:,:)
40    INTEGER :: ind
41   
42    R0=radius*0.5
43    rt=radius*0.5
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  CHARACTER(len=255) :: ncar_adv_shape
81
82  u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0
83 
84  DO l=1, llm+1
85     pr = ap(l) + bp(l)*ncar_p0
86     zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0) 
87  ENDDO
88 
89  DO l=1, llm 
90     zrl(l) = 0.5*(zr(l) + zr(l+1)) 
91  END DO
92 
93  ncar_adv_shape='cos_bell'
94  CALL getin('ncar_adv_shape',ncar_adv_shape)
95 
96  SELECT CASE(TRIM(ncar_adv_shape))
97     !--------------------------------------------- SINGLE COSINE BELL
98  CASE('const')
99     q=1
100  CASE('cos_bell')
101     CALL cosine_bell_1(q) 
102     
103  CASE('slotted_cyl')
104     CALL slotted_cylinders(q) 
105     
106  CASE('dbl_cos_bell_q1')
107     CALL cosine_bell_2(q) 
108     
109  CASE('dbl_cos_bell_q2')
110     CALL cosine_bell_2(q) 
111     DO l=1,llm
112        q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l)
113     END DO
114     
115  CASE('complement')
116     ! tracer such that, in combination with the other tracer fields
117     ! with weight (3/10), the sum is equal to one
118     CALL cosine_bell_2(qxt1)
119     DO l = 1,llm
120        q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l) 
121     END DO
122     q = q + qxt1 
123     CALL slotted_cylinders(qxt1)
124     q = q + qxt1 
125     q = 1. - q*0.3 
126     
127  CASE('hadley')  ! hadley like meridional circulation
128     CALL hadleyq(q) 
129     
130  CASE DEFAULT
131     PRINT *, 'Bad selector for variable ncar_adv_shape : <', TRIM(ncar_adv_shape),   &
132          '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
133          '<dbl_cos_bell_q2>, <complement>, <hadley>'
134     STOP
135     
136  END SELECT
137     
138CONTAINS
139
140!======================================================================
141
142    SUBROUTINE cosine_bell_1(hx)
143    IMPLICIT NONE
144    REAL(rstd) :: hx(iim*jjm,llm)
145   
146      DO l=1,llm 
147        DO j=jj_begin-1,jj_end+1
148          DO i=ii_begin-1,ii_end+1
149            n=(j-1)*iim+i
150            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
151            CALL dist_lonlat(lon0,lat0,lon,lat,rr1)  ! GC distance from center
152            rr1 = radius*rr1
153            IF ( rr1 .LT. R0 ) then
154              hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
155            ELSE
156              hx(n,l)=0.0
157            END IF
158          END DO
159        END DO
160      END DO
161   
162    END SUBROUTINE cosine_bell_1
163
164!==============================================================================
165
166    SUBROUTINE   cosine_bell_2(hx) 
167    IMPLICIT NONE       
168    REAL(rstd) :: hx(iim*jjm,llm)
169   
170      DO l=1,llm
171        DO j=jj_begin-1,jj_end+1
172          DO i=ii_begin-1,ii_end+1
173            n=(j-1)*iim+i
174            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
175            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center
176            rr1 = radius*rr1
177            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center
178            rr2 = radius*rr2
179            dd1t1 = rr1/rt
180            dd1t1 = dd1t1*dd1t1
181            dd1t2 = (zrl(l) - zc)/zt
182            dd1t2 = dd1t2*dd1t2 
183            dd1 = dd1t1 + dd1t2 
184            dd1 = Min(1.0,dd1) 
185            dd2t1 = rr2/rt
186            dd2t1 = dd2t1*dd2t1
187            dd2 = dd2t1 + dd1t2 
188            dd2 = Min(1.0,dd2) 
189            hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
190          END DO
191        END DO
192      END DO
193
194    END SUBROUTINE cosine_bell_2
195
196!=============================================================================
197
198    SUBROUTINE slotted_cylinders(hx) 
199    IMPLICIT NONE
200    REAL(rstd) :: hx(iim*jjm,llm)
201   
202      DO l=1,llm
203        DO j=jj_begin-1,jj_end+1
204          DO i=ii_begin-1,ii_end+1
205            n=(j-1)*iim+i
206            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
207            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center
208            rr1 = radius*rr1
209            CALL dist_lonlat(lonc2,latc2,lon,lat,rr2)  ! GC distance from center
210            rr2 = radius*rr2
211            dd1t1 = rr1/rt
212            dd1t1 = dd1t1*dd1t1
213            dd1t2 = (zrl(l) - zc)/zt
214            dd1t2 = dd1t2*dd1t2 
215            dd1 = dd1t1 + dd1t2 
216            dd2t1 = rr2/rt
217            dd2t1 = dd2t1*dd2t1
218            dd2 = dd2t1 + dd1t2
219           
220            IF ( dd1 .LT. 0.5 ) Then
221              hx(n,l) = 1.0
222            ELSEIF ( dd2 .LT. 0.5 ) Then
223              hx(n,l) = 1.0
224            ELSE
225              hx(n,l) = 0.1 
226            END IF 
227   
228            IF ( zrl(l) .GT. zc ) Then
229              IF ( ABS(latc1 - lat) .LT. 0.125 ) Then
230                hx(n,l)= 0.0
231              ENDIF
232            ENDIF 
233 
234          ENDDO
235       END DO
236      END DO
237
238    END SUBROUTINE slotted_cylinders 
239
240!==============================================================================
241
242    SUBROUTINE hadleyq(hx) 
243    IMPLICIT NONE
244    REAL(rstd)::hx(iim*jjm,llm) 
245    REAL(rstd),PARAMETER:: zz1=3500.,zz2=6500.,zz0=0.5*(zz1+zz2)
246   
247      DO l=1,llm
248        IF ( ( zz1 .LT. zrl(l) ) .and. ( zrl(l) .LT. zz2 ) )  THEN
249          hx(:,l) = 0.5*(1. + cos(0.002*pi*(zrl(l)-zz0)/3.)) 
250        ELSE
251          hx(:,l) = 0.0 
252        END IF
253      END DO
254    END SUBROUTINE hadleyq 
255
256  END SUBROUTINE compute_etat0_ncar
257 
258
259END MODULE etat0_ncar_mod
Note: See TracBrowser for help on using the repository browser.