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

Last change on this file since 1048 was 1048, checked in by ymipsl, 4 years ago

Introduce modification from A. Durocher github to make held&suarez testcase working on GPU
==> missing kernel from rev. 1046
YM & AD
File size: 7.7 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 :: knewt_t(:),kfrict(:)
11!$OMP THREADPRIVATE(knewt_t,kfrict)
12  LOGICAL :: done=.FALSE.
13!$OMP THREADPRIVATE(done)
14
15  REAL(rstd) :: p0,teta0,ttp,delt_y,delt_z,eps
16!$OMP THREADPRIVATE(p0,teta0,ttp,delt_y,delt_z,eps)
17
18  REAL(rstd) :: 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    CALL allocate_field(f_ps,field_t,type_real)
35    CALL allocate_field(f_phis,field_t,type_real)
36    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,2)
37    CALL allocate_field(f_u,field_u,type_real,llm)
38    CALL allocate_field(f_Ki,field_t,type_real,llm)
39
40    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
41    CALL kinetic(f_u,f_Ki)
42
43    CALL writefield('ps',f_ps)
44    CALL writefield('theta',f_theta_rhodz)
45  END SUBROUTINE test_etat0_heldsz
46
47  SUBROUTINE init_etat0
48    p0=1e5 ! default value of initial surface pressure as in H&S paper
49    ! p0=101080 for CMIP5 aquaplanets, cf LMDZ5 ini_aqua
50    CALL getin('heldsz_p0',p0)
51  END SUBROUTINE init_etat0
52
53  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
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 :: theta_eq(:,:) 
68    REAL(rstd),POINTER :: theta(:,:) 
69
70    INTEGER :: ind
71
72    CALL Init_Teq
73    DO ind=1,ndomain
74       IF (.NOT. assigned_domain(ind)) CYCLE
75       CALL swap_dimensions(ind)
76       CALL swap_geometry(ind)
77
78       theta_eq=f_theta_eq(ind) 
79       CALL compute_Teq(lat_i,theta_eq) ! FIXME : already done by Init_Teq
80
81       ps=f_ps(ind)
82       phis=f_phis(ind)
83       u=f_u(ind)
84       ps(:)=p0
85       phis(:)=0.
86       u(:,:)=0
87
88       theta_rhodz=f_theta_rhodz(ind)
89       theta=f_theta(ind)
90       CALL compute_etat0_heldsz(theta_eq,theta)
91       CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz(:,:,1),1, ondevice=.false.)
92       IF(nqtot>0) THEN
93          q=f_q(ind)
94          q(:,:,1)=0.
95          IF(nqtot>1) q(:,:,2)=0
96          IF(nqtot>2) q(:,:,3:)=0.
97       END IF
98   
99       call update_device_field(f_theta_eq)
100       call update_device_field(f_theta)
101      !$acc enter data copyin(knewt_t(:)) async
102      !$acc enter data copyin(kfrict(:)) async
103
104    ENDDO
105  END SUBROUTINE etat0
106
107  SUBROUTINE init_Teq
108    USE abort_mod
109    USE disvert_mod, ONLY : ap,bp
110    REAL(rstd),POINTER :: theta_eq(:,:) 
111    REAL(rstd) :: zsig
112    INTEGER :: ind, l
113
114    IF(.NOT.done) THEN
115       done = .TRUE.
116       
117       CALL allocate_field(f_theta,field_t,type_real,llm)
118       CALL allocate_field(f_theta_eq,field_t,type_real,llm)
119       ALLOCATE(knewt_t(llm)); ALLOCATE( kfrict(llm)) 
120
121       k_f=1.                !friction
122       CALL getin('k_j',k_f)
123       k_f=1./(daysec*k_f)
124       k_c_s=4.  !cooling surface
125       CALL getin('k_c_s',k_c_s)
126       k_c_s=1./(daysec*k_c_s)
127       k_c_a=40. !cooling free atm
128       CALL getin('k_c_a',k_c_a)
129       k_c_a=1./(daysec*k_c_a)
130       ! Constants for Teta equilibrium profile
131       teta0=315.     ! mean Teta (S.H. 315K)
132       CALL getin('teta0',teta0)
133       ttp=200.       ! Tropopause temperature (S.H. 200K)
134       CALL getin('ttp',ttp)
135       eps=0.         ! Deviation to N-S symmetry(~0-20K)
136       CALL getin('eps',eps)
137       delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
138       CALL getin('delt_y',delt_y)
139       delt_z=10.     ! Vertical Gradient (S.H. 10K)
140       CALL getin('delt_z',delt_z)
141
142       !-----------------------------------------------------------
143       knewt_g=k_c_a 
144       DO l=1,llm
145          zsig=ap(l)/preff+bp(l)
146          knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) 
147          kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) 
148       ENDDO
149
150       DO ind=1,ndomain
151          IF (.NOT. assigned_domain(ind)) CYCLE
152          CALL swap_dimensions(ind)
153          CALL swap_geometry(ind)
154          theta_eq=f_theta_eq(ind)
155          CALL compute_Teq(lat_i,theta_eq)
156       ENDDO
157
158    ELSE
159       CALL dynamico_abort( "Init_Teq called twice" )
160    END IF
161
162  END SUBROUTINE init_Teq
163
164  SUBROUTINE compute_Teq(lat,theta_eq)
165    USE disvert_mod
166    REAL(rstd),INTENT(IN) :: lat(iim*jjm)
167    REAL(rstd),INTENT(OUT) :: theta_eq(iim*jjm,llm) 
168
169    REAL(rstd) :: zsig, ddsin, tetastrat, tetajl
170    INTEGER :: i,j,l,ij
171
172    DO l=1,llm
173       zsig=ap(l)/preff+bp(l)
174       tetastrat=ttp*zsig**(-kappa)
175       DO j=jj_begin-1,jj_end+1
176          DO i=ii_begin-1,ii_end+1
177             ij=(j-1)*iim+i
178             ddsin=SIN(lat(ij)) 
179             tetajl=teta0-delt_y*ddsin*ddsin+eps*ddsin &
180                  -delt_z*(1.-ddsin*ddsin)*log(zsig)
181             theta_eq(ij,l)=MAX(tetajl,tetastrat) 
182          ENDDO
183       ENDDO
184    ENDDO
185  END SUBROUTINE compute_Teq
186
187  SUBROUTINE compute_etat0_heldsz(theta_eq, theta)
188    USE disvert_mod
189    REAL(rstd),INTENT(IN) :: theta_eq(iim*jjm,llm) 
190    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm) 
191
192    REAL(rstd) :: r  ! random number
193    INTEGER :: i,j,l,ij
194
195    DO l=1,llm
196       DO j=jj_begin-1,jj_end+1
197          DO i=ii_begin-1,ii_end+1
198             ij=(j-1)*iim+i
199             CALL RANDOM_NUMBER(r); r = 0.0 
200             theta(ij,l)=theta_eq(ij,l)*(1.+0.0005*r)
201          ENDDO
202       ENDDO
203    ENDDO
204
205  END SUBROUTINE compute_etat0_heldsz
206
207
208  SUBROUTINE held_suarez(f_ps,f_theta_rhodz,f_u) 
209    TYPE(t_field),POINTER :: f_theta_rhodz(:)
210    TYPE(t_field),POINTER :: f_u(:)
211    TYPE(t_field),POINTER :: f_ps(:)
212    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
213    REAL(rstd),POINTER :: u(:,:)
214    REAL(rstd),POINTER :: ps(:)
215    REAL(rstd),POINTER :: theta_eq(:,:)
216    REAL(rstd),POINTER :: theta(:,:)
217    INTEGER::ind
218
219    DO ind=1,ndomain
220       IF (.NOT. assigned_domain(ind)) CYCLE
221       CALL swap_dimensions(ind)
222       CALL swap_geometry(ind)
223       theta_rhodz=f_theta_rhodz(ind)
224       u=f_u(ind)
225       ps=f_ps(ind) 
226       theta_eq=f_theta_eq(ind) 
227       theta=f_theta(ind) 
228       CALL compute_heldsz(ps,theta_eq,lat_i, theta_rhodz(:,:,1),u, theta) 
229    ENDDO
230  END SUBROUTINE held_suarez
231
232  SUBROUTINE compute_heldsz(ps,theta_eq,lat, theta_rhodz,u, theta) 
233    USE theta2theta_rhodz_mod
234    USE omp_para
235    REAL(rstd),INTENT(IN)    :: ps(iim*jjm) 
236    REAL(rstd),INTENT(IN)    :: theta_eq(iim*jjm,llm) 
237    REAL(rstd),INTENT(IN)    :: lat(iim*jjm) 
238    REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm)
239    REAL(rstd),INTENT(INOUT) :: u(3*iim*jjm,llm)
240    REAL(rstd),INTENT(OUT)   :: theta(iim*jjm,llm) 
241
242    INTEGER :: i,j,l,ij
243
244    CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1, ondevice=.TRUE.)
245
246    !$acc parallel loop collapse(3) default(present) async
247    DO l=ll_begin,ll_end
248       DO j=jj_begin-1,jj_end+1
249          DO i=ii_begin-1,ii_end+1
250             ij=(j-1)*iim+i
251             theta(ij,l)=theta(ij,l) - itau_physics*dt*(theta(ij,l)-theta_eq(ij,l))* &
252                  (knewt_g+knewt_t(l)*COS(lat(ij))**4 )
253          ENDDO
254       ENDDO
255    ENDDO
256    CALL compute_theta2theta_rhodz(ps,theta,theta_rhodz,1, ondevice=.true.)
257
258    !$acc kernels default(present) async
259    DO l=ll_begin,ll_end
260       u(:,l)=u(:,l)*(1.-itau_physics*dt*kfrict(l))
261    END DO
262    !$acc end kernels
263
264  END SUBROUTINE compute_heldsz
265
266END MODULE etat0_heldsz_mod
Note: See TracBrowser for help on using the repository browser.