source: codes/icosagcm/trunk/src/initial/etat0_heldsz.f90 @ 548

Last change on this file since 548 was 548, checked in by dubos, 7 years ago

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