source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/etat0_heldsz.f90 @ 316

Last change on this file since 316 was 221, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM
File size: 7.8 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,SAVE :: knewt_t(:),kfrict(:)
11!$OMP THREADPRIVATE(knewt_t,kfrict)
12  LOGICAL, SAVE :: done=.FALSE.
13!$OMP THREADPRIVATE(done)
14
15  REAL(rstd),SAVE :: teta0,ttp,delt_y,delt_z,eps
16!$OMP THREADPRIVATE(teta0,ttp,delt_y,delt_z,eps)
17
18  REAL(rstd),SAVE :: knewt_g, k_f,k_c_a,k_c_s
19!$OMP THREADPRIVATE(knewt_g, k_f,k_c_a,k_c_s)
20
21  PUBLIC :: etat0, held_suarez
22 
23CONTAINS
24
25  SUBROUTINE test_etat0_heldsz
26    USE icosa
27    USE kinetic_mod
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    TYPE(t_field),POINTER :: f_Ki(:)
35
36    REAL(rstd),POINTER :: Ki(:,:)
37    INTEGER :: ind
38
39    CALL allocate_field(f_ps,field_t,type_real)
40    CALL allocate_field(f_phis,field_t,type_real)
41    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
42    CALL allocate_field(f_u,field_u,type_real,llm)
43    CALL allocate_field(f_Ki,field_t,type_real,llm)
44
45    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
46    CALL kinetic(f_u,f_Ki)
47
48    CALL writefield('ps',f_ps)
49    CALL writefield('theta',f_theta_rhodz)
50  END SUBROUTINE test_etat0_heldsz
51
52  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
53    USE icosa
54    USE theta2theta_rhodz_mod
55    IMPLICIT NONE
56    TYPE(t_field),POINTER :: f_ps(:)
57    TYPE(t_field),POINTER :: f_phis(:)
58    TYPE(t_field),POINTER :: f_theta_rhodz(:)
59    TYPE(t_field),POINTER :: f_u(:)
60    TYPE(t_field),POINTER :: f_q(:)
61
62    REAL(rstd),POINTER :: ps(:)
63    REAL(rstd),POINTER :: phis(:)
64    REAL(rstd),POINTER :: theta_rhodz(:,:)
65    REAL(rstd),POINTER :: u(:,:)
66    REAL(rstd),POINTER :: q(:,:,:)
67    REAL(rstd),POINTER :: clat(:) 
68    REAL(rstd),POINTER :: theta_eq(:,:) 
69    REAL(rstd),POINTER :: theta(:,:) 
70
71    INTEGER :: ind
72
73    CALL Init_Teq
74    DO ind=1,ndomain
75       IF (.NOT. assigned_domain(ind)) CYCLE
76       CALL swap_dimensions(ind)
77       CALL swap_geometry(ind)
78
79       theta_eq=f_theta_eq(ind) 
80       clat=f_clat(ind) 
81       CALL compute_Teq(clat,theta_eq) ! FIXME : already done by Init_Teq
82
83       ps=f_ps(ind)
84       phis=f_phis(ind)
85       u=f_u(ind)
86       ps(:)=1e5
87       phis(:)=0.
88       u(:,:)=0
89
90       theta_rhodz=f_theta_rhodz(ind)
91       theta=f_theta(ind)
92       CALL compute_etat0_heldsz(theta_eq,theta)
93       CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
94       q=f_q(ind)
95       q(:,:,:)=1e2
96    ENDDO
97  END SUBROUTINE etat0
98
99  SUBROUTINE init_Teq
100    USE icosa
101    USE disvert_mod, ONLY : ap,bp
102    IMPLICIT NONE
103    REAL(rstd),POINTER :: clat(:) 
104    REAL(rstd),POINTER :: theta_eq(:,:) 
105    REAL(rstd) :: zsig
106    INTEGER :: ind, l
107
108    IF(.NOT.done) THEN
109       done = .TRUE.
110       
111       CALL allocate_field(f_theta,field_t,type_real,llm)
112       CALL allocate_field(f_theta_eq,field_t,type_real,llm)
113       CALL allocate_field(f_clat,field_t,type_real)
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          clat=f_clat(ind)
150          theta_eq=f_theta_eq(ind)
151          CALL compute_Teq(clat,theta_eq)
152       ENDDO
153
154    ELSE
155       PRINT *, 'Init_Teq called twice'
156       CALL ABORT
157    END IF
158
159  END SUBROUTINE init_Teq
160
161  SUBROUTINE compute_Teq(clat,theta_eq)
162    USE icosa
163    USE disvert_mod
164    IMPLICIT NONE 
165    REAL(rstd),INTENT(OUT) :: clat(iim*jjm)
166    REAL(rstd),INTENT(OUT) :: theta_eq(iim*jjm,llm) 
167
168    REAL(rstd) :: lon, lat, r, zsig, ddsin, tetastrat, tetajl
169    REAL(rstd) :: slat(iim*jjm) 
170    INTEGER :: i,j,l,ij
171
172    DO j=jj_begin-1,jj_end+1
173       DO i=ii_begin-1,ii_end+1
174          ij=(j-1)*iim+i
175          CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
176          clat(ij)=cos(lat) 
177          slat(ij)=sin(lat) 
178       ENDDO
179    ENDDO
180
181    DO l=1,llm
182       zsig=ap(l)/preff+bp(l)
183       tetastrat=ttp*zsig**(-kappa)
184       DO j=jj_begin-1,jj_end+1
185          DO i=ii_begin-1,ii_end+1
186             ij=(j-1)*iim+i
187             ddsin=slat(ij) 
188             tetajl=teta0-delt_y*ddsin*ddsin+eps*ddsin &
189                  -delt_z*(1.-ddsin*ddsin)*log(zsig)
190             theta_eq(ij,l)=MAX(tetajl,tetastrat) 
191          ENDDO
192       ENDDO
193    ENDDO
194  END SUBROUTINE compute_Teq
195
196  SUBROUTINE compute_etat0_heldsz(theta_eq, theta)
197    USE icosa
198    USE disvert_mod
199    USE pression_mod
200    USE exner_mod
201    USE geopotential_mod
202    USE theta2theta_rhodz_mod
203    IMPLICIT NONE 
204    REAL(rstd),INTENT(IN) :: theta_eq(iim*jjm,llm) 
205    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
206
207    REAL(rstd) :: r  ! random number
208    INTEGER :: i,j,l,ij
209
210    DO l=1,llm
211       DO j=jj_begin-1,jj_end+1
212          DO i=ii_begin-1,ii_end+1
213             ij=(j-1)*iim+i
214             CALL RANDOM_NUMBER(r); r = 0.0 
215             theta(ij,l)=theta_eq(ij,l)*(1.+0.0005*r)
216          ENDDO
217       ENDDO
218    ENDDO
219
220  END SUBROUTINE compute_etat0_heldsz
221
222
223  SUBROUTINE held_suarez(f_ps,f_theta_rhodz,f_u) 
224    USE icosa
225    IMPLICIT NONE
226    TYPE(t_field),POINTER :: f_theta_rhodz(:)
227    TYPE(t_field),POINTER :: f_u(:)
228    TYPE(t_field),POINTER :: f_ps(:)
229    REAL(rstd),POINTER :: theta_rhodz(:,:)
230    REAL(rstd),POINTER :: u(:,:)
231    REAL(rstd),POINTER :: ps(:)
232    REAL(rstd),POINTER :: theta_eq(:,:)
233    REAL(rstd),POINTER :: theta(:,:)
234    REAL(rstd),POINTER :: clat(:)
235    INTEGER::ind
236
237    DO ind=1,ndomain
238       IF (.NOT. assigned_domain(ind)) CYCLE
239       CALL swap_dimensions(ind)
240       CALL swap_geometry(ind)
241       theta_rhodz=f_theta_rhodz(ind)
242       u=f_u(ind)
243       ps=f_ps(ind) 
244       theta_eq=f_theta_eq(ind) 
245       theta=f_theta(ind) 
246       clat=f_clat(ind) 
247       CALL compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta) 
248    ENDDO
249  END SUBROUTINE held_suarez
250
251  SUBROUTINE compute_heldsz(ps,theta_eq,clat, theta_rhodz,u, theta) 
252    USE icosa
253    USE theta2theta_rhodz_mod
254    IMPLICIT NONE
255    REAL(rstd),INTENT(IN)    :: ps(iim*jjm) 
256    REAL(rstd),INTENT(IN)    :: theta_eq(iim*jjm,llm) 
257    REAL(rstd),INTENT(IN)    :: clat(iim*jjm) 
258    REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm)
259    REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm)
260    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm) 
261
262    INTEGER :: i,j,l,ij
263
264    CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
265    DO l=1,llm
266       DO j=jj_begin-1,jj_end+1
267          DO i=ii_begin-1,ii_end+1
268             ij=(j-1)*iim+i
269             theta(ij,l)=theta(ij,l) - dt*(theta(ij,l)-theta_eq(ij,l))* &
270                  (knewt_g+knewt_t(l)*clat(ij)**4 )
271          ENDDO
272       ENDDO
273    ENDDO
274    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1)
275
276    Do l=1,llm
277       u(:,l)=u(:,l)*(1.-dt*kfrict(l))
278    END DO
279
280  END SUBROUTINE compute_heldsz
281
282END MODULE etat0_heldsz_mod
Note: See TracBrowser for help on using the repository browser.