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

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

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