source: codes/icosagcm/trunk/src/etat0_heldsz.f90 @ 170

Last change on this file since 170 was 170, checked in by dubos, 11 years ago

Activated call to physics - Held & Suarez test case seems to work now
File size: 7.5 KB
Line 
1MODULE etat0_heldsz_mod
2  USE icosa
3  IMPLICIT NONE
4  PRIVATE
5
6  TYPE(t_field),POINTER :: f_theta_eq(:)
7  TYPE(t_field),POINTER :: f_theta(:)
8  TYPE(t_field),POINTER :: f_clat(:) ! FIXME, duplication
9
10  REAL(rstd),ALLOCATABLE :: knewt_t(:),kfrict(:)
11
12  LOGICAL, SAVE :: done=.FALSE.
13
14  REAL(rstd) :: teta0,ttp,delt_y,delt_z,eps
15  REAL(rstd) :: knewt_g, k_f,k_c_a,k_c_s
16
17  PUBLIC :: etat0, held_suarez
18 
19CONTAINS
20
21  SUBROUTINE test_etat0_heldsz
22    USE icosa
23    USE kinetic_mod
24    IMPLICIT NONE
25    TYPE(t_field),POINTER :: f_ps(:)
26    TYPE(t_field),POINTER :: f_phis(:)
27    TYPE(t_field),POINTER :: f_theta_rhodz(:)
28    TYPE(t_field),POINTER :: f_u(:)
29    TYPE(t_field),POINTER :: f_q(:)
30    TYPE(t_field),POINTER :: f_Ki(:)
31
32    REAL(rstd),POINTER :: Ki(:,:)
33    INTEGER :: ind
34
35    CALL allocate_field(f_ps,field_t,type_real)
36    CALL allocate_field(f_phis,field_t,type_real)
37    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
38    CALL allocate_field(f_u,field_u,type_real,llm)
39    CALL allocate_field(f_Ki,field_t,type_real,llm)
40
41    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
42    CALL kinetic(f_u,f_Ki)
43
44    CALL writefield('ps',f_ps)
45    CALL writefield('theta',f_theta_rhodz)
46  END SUBROUTINE test_etat0_heldsz
47
48  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
49    USE icosa
50    USE theta2theta_rhodz_mod
51    IMPLICIT NONE
52    TYPE(t_field),POINTER :: f_ps(:)
53    TYPE(t_field),POINTER :: f_phis(:)
54    TYPE(t_field),POINTER :: f_theta_rhodz(:)
55    TYPE(t_field),POINTER :: f_u(:)
56    TYPE(t_field),POINTER :: f_q(:)
57
58    REAL(rstd),POINTER :: ps(:)
59    REAL(rstd),POINTER :: phis(:)
60    REAL(rstd),POINTER :: theta_rhodz(:,:)
61    REAL(rstd),POINTER :: u(:,:)
62    REAL(rstd),POINTER :: q(:,:,:)
63    REAL(rstd),POINTER :: clat(:) 
64    REAL(rstd),POINTER :: theta_eq(:,:) 
65    REAL(rstd),POINTER :: theta(:,:) 
66
67    INTEGER :: ind
68
69    CALL Init_Teq
70    DO ind=1,ndomain
71       CALL swap_dimensions(ind)
72       CALL swap_geometry(ind)
73
74       theta_eq=f_theta_eq(ind) 
75       clat=f_clat(ind) 
76       CALL compute_Teq(clat,theta_eq) ! FIXME : already done by Init_Teq
77
78       ps=f_ps(ind)
79       phis=f_phis(ind)
80       u=f_u(ind)
81       ps(:)=1e5
82       phis(:)=0.
83       u(:,:)=0
84
85       theta_rhodz=f_theta_rhodz(ind)
86       theta=f_theta(ind)
87       CALL compute_etat0_heldsz(theta_eq,theta)
88       CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
89       q=f_q(ind)
90       q(:,:,:)=1e2
91    ENDDO
92  END SUBROUTINE etat0
93
94  SUBROUTINE init_Teq
95    USE icosa
96    USE disvert_mod, ONLY : ap,bp
97    IMPLICIT NONE
98    REAL(rstd),POINTER :: clat(:) 
99    REAL(rstd),POINTER :: theta_eq(:,:) 
100    REAL(rstd) :: zsig
101    INTEGER :: ind, l
102
103    IF(.NOT.done) THEN
104       done = .TRUE.
105       
106       CALL allocate_field(f_theta,field_t,type_real,llm)
107       CALL allocate_field(f_theta_eq,field_t,type_real,llm)
108       CALL allocate_field(f_clat,field_t,type_real)
109       ALLOCATE(knewt_t(llm)); ALLOCATE( kfrict(llm)) 
110
111       k_f=1.                !friction
112       CALL getin('k_j',k_f)
113       k_f=1./(daysec*k_f)
114       k_c_s=4.  !cooling surface
115       CALL getin('k_c_s',k_c_s)
116       k_c_s=1./(daysec*k_c_s)
117       k_c_a=40. !cooling free atm
118       CALL getin('k_c_a',k_c_a)
119       k_c_a=1./(daysec*k_c_a)
120       ! Constants for Teta equilibrium profile
121       teta0=315.     ! mean Teta (S.H. 315K)
122       CALL getin('teta0',teta0)
123       ttp=200.       ! Tropopause temperature (S.H. 200K)
124       CALL getin('ttp',ttp)
125       eps=0.         ! Deviation to N-S symmetry(~0-20K)
126       CALL getin('eps',eps)
127       delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
128       CALL getin('delt_y',delt_y)
129       delt_z=10.     ! Vertical Gradient (S.H. 10K)
130       CALL getin('delt_z',delt_z)
131
132       !-----------------------------------------------------------
133       knewt_g=k_c_a 
134       DO l=1,llm
135          zsig=ap(l)/preff+bp(l)
136          knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 
137          kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 
138       ENDDO
139
140       DO ind=1,ndomain
141          CALL swap_dimensions(ind)
142          CALL swap_geometry(ind)
143          clat=f_clat(ind)
144          theta_eq=f_theta_eq(ind)
145          CALL compute_Teq(clat,theta_eq)
146       ENDDO
147
148    ELSE
149       PRINT *, 'Init_Teq called twice'
150       CALL ABORT
151    END IF
152
153  END SUBROUTINE init_Teq
154
155  SUBROUTINE compute_Teq(clat,theta_eq)
156    USE icosa
157    USE disvert_mod
158    IMPLICIT NONE 
159    REAL(rstd),INTENT(OUT) :: clat(iim*jjm)
160    REAL(rstd),INTENT(OUT) :: theta_eq(iim*jjm,llm) 
161
162    REAL(rstd) :: lon, lat, r, zsig, ddsin, tetastrat, tetajl
163    REAL(rstd) :: slat(iim*jjm) 
164    INTEGER :: i,j,l,ij
165
166    DO j=jj_begin-1,jj_end+1
167       DO i=ii_begin-1,ii_end+1
168          ij=(j-1)*iim+i
169          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
170          clat(ij)=cos(lat) 
171          slat(ij)=sin(lat) 
172       ENDDO
173    ENDDO
174
175    DO l=1,llm
176       zsig=ap(l)/preff+bp(l)
177       tetastrat=ttp*zsig**(-kappa)
178       DO j=jj_begin-1,jj_end+1
179          DO i=ii_begin-1,ii_end+1
180             ij=(j-1)*iim+i
181             ddsin=slat(ij) 
182             tetajl=teta0-delt_y*ddsin*ddsin+eps*ddsin &
183                  -delt_z*(1.-ddsin*ddsin)*log(zsig)
184             theta_eq(ij,l)=MAX(tetajl,tetastrat) 
185          ENDDO
186       ENDDO
187    ENDDO
188  END SUBROUTINE compute_Teq
189
190  SUBROUTINE compute_etat0_heldsz(theta_eq, theta)
191    USE icosa
192    USE disvert_mod
193    USE pression_mod
194    USE exner_mod
195    USE geopotential_mod
196    USE theta2theta_rhodz_mod
197    IMPLICIT NONE 
198    REAL(rstd),INTENT(IN) :: theta_eq(iim*jjm,llm) 
199    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
200
201    REAL(rstd) :: r  ! random number
202    INTEGER :: i,j,l,ij
203
204    DO l=1,llm
205       DO j=jj_begin-1,jj_end+1
206          DO i=ii_begin-1,ii_end+1
207             ij=(j-1)*iim+i
208             CALL RANDOM_NUMBER(r); r = 0.0 
209             theta(ij,l)=theta_eq(ij,l)*(1.+0.0005*r)
210          ENDDO
211       ENDDO
212    ENDDO
213
214  END SUBROUTINE compute_etat0_heldsz
215
216
217  SUBROUTINE held_suarez(f_ps,f_theta_rhodz,f_u) 
218    USE icosa
219    IMPLICIT NONE
220    TYPE(t_field),POINTER :: f_theta_rhodz(:)
221    TYPE(t_field),POINTER :: f_u(:)
222    TYPE(t_field),POINTER :: f_ps(:)
223    REAL(rstd),POINTER :: theta_rhodz(:,:)
224    REAL(rstd),POINTER :: u(:,:)
225    REAL(rstd),POINTER :: ps(:)
226    REAL(rstd),POINTER :: theta_eq(:,:)
227    REAL(rstd),POINTER :: theta(:,:)
228    REAL(rstd),POINTER :: clat(:)
229    INTEGER::ind
230
231    DO ind=1,ndomain
232       CALL swap_dimensions(ind)
233       CALL swap_geometry(ind)
234       theta_rhodz=f_theta_rhodz(ind)
235       u=f_u(ind)
236       ps=f_ps(ind) 
237       theta_eq=f_theta_eq(ind) 
238       theta=f_theta(ind) 
239       clat=f_clat(ind) 
240       CALL compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta) 
241    ENDDO
242  END SUBROUTINE held_suarez
243
244  SUBROUTINE compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta) 
245    USE icosa
246    USE theta2theta_rhodz_mod
247    IMPLICIT NONE
248    REAL(rstd),INTENT(IN)    :: ps(iim*jjm) 
249    REAL(rstd),INTENT(IN)    :: theta_eq(iim*jjm,llm) 
250    REAL(rstd),INTENT(IN)    :: clat(iim*jjm) 
251    REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm)
252    REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm)
253    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm) 
254
255    INTEGER :: i,j,l,ij
256
257    CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
258    DO l=1,llm
259       DO j=jj_begin-1,jj_end+1
260          DO i=ii_begin-1,ii_end+1
261             ij=(j-1)*iim+i
262             theta(ij,l)=theta(ij,l) - dt*(theta(ij,l)-theta_eq(ij,l))* &
263                  (knewt_g+knewt_t(l)*clat(ij)**4 )
264          ENDDO
265       ENDDO
266    ENDDO
267    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
268
269    Do l=1,llm
270       u(:,l)=u(:,l)*(1.-dt*kfrict(l))
271    END DO
272
273  END SUBROUTINE compute_heldsz
274
275END MODULE etat0_heldsz_mod
Note: See TracBrowser for help on using the repository browser.