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

Last change on this file since 344 was 344, checked in by dubos, 9 years ago

Cleanup DCMIP 1&2
File size: 7.4 KB
Line 
1MODULE etat0_dcmip1_mod
2  USE icosa
3  IMPLICIT NONE         
4  PRIVATE
5
6  INTEGER, PARAMETER :: const=1, cos_bell=2, slotted_cyl=3, &
7       dbl_cos_bell_q1=4, dbl_cos_bell_q2=5, complement=6, hadley=7, &
8       dcmip11=-1
9  INTEGER, SAVE :: shape
10!$OMP THREADPRIVATE(shape)
11
12  REAL(rstd), SAVE  :: h0=1.
13!$OMP THREADPRIVATE(h0)
14  REAL(rstd), SAVE  :: lon0=3*pi/2
15!$OMP THREADPRIVATE(lon0)
16  REAL(rstd), SAVE  :: lat0=0.0
17!$OMP THREADPRIVATE(lat0)
18  REAL(rstd), SAVE  :: alpha=0.0
19!$OMP THREADPRIVATE(alpha)
20  REAL(rstd), SAVE  :: R0 
21!$OMP THREADPRIVATE(R0)
22  REAL(rstd), SAVE  :: lat1=0.
23!$OMP THREADPRIVATE(lat1)
24  REAL(rstd), SAVE  :: lat2=0.
25!$OMP THREADPRIVATE(lat2)
26  REAL(rstd), SAVE  :: lon1=pi/6
27!$OMP THREADPRIVATE(lon1)
28  REAL(rstd), SAVE  :: lon2=-pi/6
29!$OMP THREADPRIVATE(lon2)
30  REAL(rstd), SAVE  :: latc1=0.
31!$OMP THREADPRIVATE(latc1)
32  REAL(rstd), SAVE  :: latc2=0.
33!$OMP THREADPRIVATE(latc2)
34  REAL(rstd), SAVE  :: lonc1=5*pi/6
35!$OMP THREADPRIVATE(lonc1)
36  REAL(rstd), SAVE  :: lonc2=7*pi/6
37!$OMP THREADPRIVATE(lonc2)
38  REAL(rstd), SAVE  :: zt=1000.0
39!$OMP THREADPRIVATE(zt)
40  REAL(rstd), SAVE  :: rt
41!$OMP THREADPRIVATE(rt)
42  REAL(rstd), SAVE  :: zc=5000.0
43!$OMP THREADPRIVATE(zc)
44
45  PUBLIC getin_etat0, compute_etat0
46 
47CONTAINS
48
49  SUBROUTINE getin_etat0
50    CHARACTER(len=255) :: dcmip1_adv_shape
51    R0=radius*0.5
52    rt=radius*0.5
53    dcmip1_adv_shape='cos_bell'
54    CALL getin('dcmip1_shape',dcmip1_adv_shape)
55    SELECT CASE(TRIM(dcmip1_adv_shape))
56    CASE('const')
57       shape=const
58    CASE('cos_bell')
59       shape=cos_bell
60    CASE('slotted_cyl')
61       shape=slotted_cyl
62    CASE('dbl_cos_bell_q1')
63       shape=dbl_cos_bell_q1
64    CASE('dbl_cos_bell_q2')
65       shape=dbl_cos_bell_q2
66    CASE('complement')
67       shape=complement
68    CASE('hadley')  ! hadley like meridional circulation
69       shape=hadley
70    CASE('dcmip11')
71       IF(nqtot<5) THEN
72          PRINT *,'Error : etat0_dcmip=dcmip11 and nqtot = ',nqtot,' .'
73          PRINT *,'nqtot must be equal to 5 when etat0_dcmip=dcmip11'
74          STOP
75       END IF
76       shape=dcmip11
77    CASE DEFAULT
78       PRINT *, 'Bad selector for variable dcmip1_adv_shape : <', TRIM(dcmip1_adv_shape),   &
79            '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
80            '<dbl_cos_bell_q2>, <complement>, <hadley>'
81       STOP
82    END SELECT
83  END SUBROUTINE getin_etat0
84
85  SUBROUTINE compute_etat0(ngrid,lon,lat, phis,ps,temp,ulon,ulat,q)
86    INTEGER, INTENT(IN) :: ngrid
87    REAL(rstd),INTENT(IN) :: lon(ngrid),lat(ngrid)
88    REAL(rstd),INTENT(OUT) :: ps(ngrid),phis(ngrid)
89    REAL(rstd),INTENT(OUT) :: temp(ngrid,llm),ulon(ngrid,llm),ulat(ngrid,llm)
90    REAL(rstd),INTENT(OUT) :: q(ngrid,llm,nqtot)
91    ps = ncar_p0
92    phis = 0. ; temp = 0. ; 
93    ulon = 0. ; ulat=0.
94    SELECT CASE(shape)
95    CASE(dcmip11)
96       CALL compute_etat0_ncar(4,ngrid,lon,lat,q(:,:,1))
97       CALL compute_etat0_ncar(5,ngrid,lon,lat,q(:,:,2))
98       CALL compute_etat0_ncar(3,ngrid,lon,lat,q(:,:,3))
99       CALL compute_etat0_ncar(6,ngrid,lon,lat,q(:,:,4))
100       CALL compute_etat0_ncar(1,ngrid,lon,lat,q(:,:,5))
101    CASE DEFAULT
102       CALL compute_etat0_ncar(shape,ngrid,lon,lat,q(:,:,1))
103    END SELECT
104  END SUBROUTINE compute_etat0
105
106  SUBROUTINE compute_etat0_ncar(icase,ngrid,lon,lat, q)
107  USE icosa
108  USE disvert_mod
109  IMPLICIT NONE 
110  INTEGER, INTENT(IN) :: icase, ngrid
111  REAL(rstd),INTENT(IN) :: lon(ngrid),lat(ngrid) 
112  REAL(rstd),INTENT(OUT) :: q(ngrid,llm) 
113  REAL(rstd) :: zr(llm+1), zrl(llm), qxt1(ngrid,llm)
114  REAL(rstd) :: pr
115  !  REAL(rstd) :: lon, lat
116  INTEGER :: n,l
117 
118  DO l=1, llm+1
119     pr = ap(l) + bp(l)*ncar_p0
120     zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0) 
121  ENDDO
122 
123  DO l=1, llm 
124     zrl(l) = 0.5*(zr(l) + zr(l+1)) 
125  END DO
126 
127  SELECT CASE(icase)
128  CASE(1)
129     q=1
130  CASE(2)
131     CALL cosine_bell_1(q)     
132  CASE(3)
133     CALL slotted_cylinders(q) 
134  CASE(4)
135     CALL cosine_bell_2(q)     
136  CASE(5)
137     CALL cosine_bell_2(q) 
138     DO l=1,llm
139        q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l)
140     END DO     
141  CASE(6)
142     ! tracer such that, in combination with the other tracer fields
143     ! with weight (3/10), the sum is equal to one
144     CALL cosine_bell_2(qxt1)
145     DO l = 1,llm
146        q(:,l) = 0.9 - 0.8*qxt1(:,l)*qxt1(:,l) 
147     END DO
148     q = q + qxt1 
149     CALL slotted_cylinders(qxt1)
150     q = q + qxt1 
151     q = 1. - q*0.3
152  CASE(7)  ! hadley like meridional circulation
153     CALL hadleyq(q) 
154  END SELECT
155     
156  CONTAINS
157!======================================================================
158
159    SUBROUTINE cosine_bell_1(hx)
160    IMPLICIT NONE
161    REAL(rstd) :: hx(ngrid,llm)
162    REAL(rstd) :: rr1,rr2   
163    INTEGER :: n,l
164    DO l=1,llm 
165       DO n=1,ngrid
166          CALL dist_lonlat(lon0,lat0,lon(n),lat(n),rr1)  ! GC distance from center
167          rr1 = radius*rr1
168          IF ( rr1 .LT. R0 ) then
169             hx(n,l)= 0.5*h0*(1+cos(pi*rr1/R0)) 
170          ELSE
171             hx(n,l)=0.0
172          END IF
173       END DO
174    END DO
175  END SUBROUTINE cosine_bell_1
176
177!==============================================================================
178 
179  SUBROUTINE cosine_bell_2(hx) 
180    REAL(rstd) :: hx(ngrid,llm)
181    REAL(rstd) :: rr1,rr2,dd1,dd2,dd1t1,dd1t2,dd2t1
182    INTEGER :: n,l
183    DO l=1,llm
184       DO n=1,ngrid
185          CALL dist_lonlat(lonc1,latc1,lon(n),lat(n),rr1)  ! GC distance from center
186          rr1 = radius*rr1
187          CALL dist_lonlat(lonc2,latc2,lon(n),lat(n),rr2)  ! GC distance from center
188          rr2 = radius*rr2
189          dd1t1 = rr1/rt
190          dd1t1 = dd1t1*dd1t1
191          dd1t2 = (zrl(l) - zc)/zt
192          dd1t2 = dd1t2*dd1t2 
193          dd1 = dd1t1 + dd1t2 
194          dd1 = Min(1.0,dd1) 
195          dd2t1 = rr2/rt
196          dd2t1 = dd2t1*dd2t1
197          dd2 = dd2t1 + dd1t2 
198          dd2 = Min(1.0,dd2) 
199          hx(n,l)= 0.5*(1. + cos(pi*dd1))+0.5*(1.+cos(pi*dd2))
200       END DO
201    END DO 
202  END SUBROUTINE cosine_bell_2
203 
204  !=============================================================================
205
206  SUBROUTINE slotted_cylinders(hx) 
207    REAL(rstd) :: hx(ngrid,llm)
208    REAL(rstd) :: rr1,rr2,dd1,dd2,dd1t1,dd1t2,dd2t1
209    INTEGER :: n,l
210    DO l=1,llm
211       DO n=1,ngrid
212          CALL dist_lonlat(lonc1,latc1,lon(n),lat(n),rr1)  ! GC distance from center
213          rr1 = radius*rr1
214          CALL dist_lonlat(lonc2,latc2,lon(n),lat(n),rr2)  ! GC distance from center
215          rr2 = radius*rr2
216          dd1t1 = rr1/rt
217          dd1t1 = dd1t1*dd1t1
218          dd1t2 = (zrl(l) - zc)/zt
219          dd1t2 = dd1t2*dd1t2 
220          dd1 = dd1t1 + dd1t2 
221          dd2t1 = rr2/rt
222          dd2t1 = dd2t1*dd2t1
223          dd2 = dd2t1 + dd1t2
224          IF ( dd1 .LT. 0.5 ) Then
225             hx(n,l) = 1.0
226          ELSEIF ( dd2 .LT. 0.5 ) Then
227             hx(n,l) = 1.0
228          ELSE
229             hx(n,l) = 0.1 
230          END IF
231          IF ( zrl(l) .GT. zc ) Then
232             IF ( ABS(latc1 - lat_i(n)) .LT. 0.125 ) Then
233                hx(n,l)= 0.1
234             ENDIF
235          ENDIF
236       END DO
237    END DO
238  END SUBROUTINE slotted_cylinders
239   
240!==============================================================================
241
242    SUBROUTINE hadleyq(hx) 
243      REAL(rstd)::hx(ngrid,llm) 
244      REAL(rstd),PARAMETER:: zz1=2000.,zz2=5000.,zz0=0.5*(zz1+zz2)
245      INTEGER :: n,l
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(2*pi*(zrl(l)-zz0)/(zz2-zz1))) 
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
258END MODULE etat0_dcmip1_mod
Note: See TracBrowser for help on using the repository browser.