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

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

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

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.