Context Navigation

source: codes/icosagcm/trunk/src/caldyn_gcm.f90 @ 409

Last change on this file since 409 was 406, checked in by dubos, 8 years ago
Introduced fake_moist mode
File size: 9.9 KB

Line
1	MODULE caldyn_gcm_mod
2	USE icosa
3	USE transfert_mod
4	USE caldyn_kernels_base_mod
5	USE caldyn_kernels_mod
6	IMPLICIT NONE
7	PRIVATE
8
9	PUBLIC init_caldyn, caldyn_BC, caldyn
10
11	CONTAINS
12
13	SUBROUTINE init_caldyn
14	USE icosa
15	USE observable_mod
16	USE mpipara
17	USE omp_para
18	IMPLICIT NONE
19	CHARACTER(len=255) :: def
20	INTEGER :: ind
21	REAL(rstd),POINTER :: planetvel(:)
22
23	def='energy'
24	CALL getin('caldyn_conserv',def)
25	SELECT CASE(TRIM(def))
26	CASE('energy')
27	caldyn_conserv=energy
28	CASE('enstrophy')
29	caldyn_conserv=enstrophy
30	CASE DEFAULT
31	IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
32	TRIM(def),'> options are <energy>, <enstrophy>'
33	STOP
34	END SELECT
35	IF (is_master) PRINT *, 'caldyn_conserv=',def
36
37	def='theta'
38	CALL getin('caldyn_thermo',def)
39	SELECT CASE(TRIM(def))
40	CASE('theta')
41	caldyn_thermo=thermo_theta
42	physics_thermo=thermo_dry
43	CASE('entropy')
44	caldyn_thermo=thermo_entropy
45	physics_thermo=thermo_dry
46	CASE('theta_fake_moist')
47	caldyn_thermo=thermo_theta
48	physics_thermo=thermo_fake_moist
49	CASE('entropy_fake_moist')
50	caldyn_thermo=thermo_entropy
51	physics_thermo=thermo_fake_moist
52	CASE('moist')
53	caldyn_thermo=thermo_moist
54	physics_thermo=thermo_moist
55	CASE DEFAULT
56	IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_thermo : <', &
57	TRIM(def),'> options are <theta>, <entropy>'
58	STOP
59	END SELECT
60
61	CALL allocate_caldyn
62
63	DO ind=1,ndomain
64	IF (.NOT. assigned_domain(ind)) CYCLE
65	CALL swap_dimensions(ind)
66	CALL swap_geometry(ind)
67	planetvel=f_planetvel(ind)
68	CALL compute_planetvel(planetvel)
69	END DO
70
71	END SUBROUTINE init_caldyn
72
73	SUBROUTINE allocate_caldyn
74	USE icosa
75	IMPLICIT NONE
76
77	CALL allocate_field(f_out_u,field_u,type_real,llm)
78	CALL allocate_field(f_qu,field_u,type_real,llm)
79	CALL allocate_field(f_qv,field_z,type_real,llm)
80	CALL allocate_field(f_pk, field_t,type_real,llm, name='pk')
81	CALL allocate_field(f_wwuu, field_u,type_real,llm+1,name='wwuu')
82	CALL allocate_field(f_planetvel, field_u,type_real, name='planetvel') ! planetary velocity at r=a
83
84	END SUBROUTINE allocate_caldyn
85
86	SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux)
87	USE icosa
88	USE mpipara
89	USE omp_para
90	TYPE(t_field),POINTER :: f_phis(:)
91	TYPE(t_field),POINTER :: f_geopot(:)
92	TYPE(t_field),POINTER :: f_wflux(:)
93	REAL(rstd),POINTER :: phis(:)
94	REAL(rstd),POINTER :: wflux(:,:)
95	REAL(rstd),POINTER :: geopot(:,:)
96	REAL(rstd),POINTER :: wwuu(:,:)
97
98	INTEGER :: ind,i,j,ij,l
99
100	IF (is_omp_first_level) THEN
101	DO ind=1,ndomain
102	IF (.NOT. assigned_domain(ind)) CYCLE
103	CALL swap_dimensions(ind)
104	CALL swap_geometry(ind)
105	geopot=f_geopot(ind)
106	phis=f_phis(ind)
107	wflux=f_wflux(ind)
108	wwuu=f_wwuu(ind)
109
110	DO ij=ij_begin_ext,ij_end_ext
111	! lower BCs : geopot=phis, wflux=0, wwuu=0
112	geopot(ij,1) = phis(ij)
113	wflux(ij,1) = 0.
114	wwuu(ij+u_right,1)=0
115	wwuu(ij+u_lup,1)=0
116	wwuu(ij+u_ldown,1)=0
117	! top BCs : wflux=0, wwuu=0
118	wflux(ij,llm+1) = 0.
119	wwuu(ij+u_right,llm+1)=0
120	wwuu(ij+u_lup,llm+1)=0
121	wwuu(ij+u_ldown,llm+1)=0
122	ENDDO
123	END DO
124	ENDIF
125
126	!$OMP BARRIER
127	END SUBROUTINE caldyn_BC
128
129	SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
130	f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
131	USE icosa
132	USE observable_mod
133	USE disvert_mod, ONLY : caldyn_eta, eta_mass
134	USE vorticity_mod
135	USE kinetic_mod
136	USE theta2theta_rhodz_mod
137	USE wind_mod
138	USE mpipara
139	USE trace
140	USE omp_para
141	USE output_field_mod
142	USE checksum_mod
143	IMPLICIT NONE
144	LOGICAL,INTENT(IN) :: write_out
145	TYPE(t_field),POINTER :: f_phis(:)
146	TYPE(t_field),POINTER :: f_ps(:)
147	TYPE(t_field),POINTER :: f_mass(:)
148	TYPE(t_field),POINTER :: f_theta_rhodz(:)
149	TYPE(t_field),POINTER :: f_u(:)
150	TYPE(t_field),POINTER :: f_q(:)
151	TYPE(t_field),POINTER :: f_geopot(:)
152	TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
153	TYPE(t_field) :: f_dps(:)
154	TYPE(t_field) :: f_dmass(:)
155	TYPE(t_field) :: f_dtheta_rhodz(:)
156	TYPE(t_field) :: f_du(:)
157
158	REAL(rstd),POINTER :: ps(:), dps(:)
159	REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
160	REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
161	REAL(rstd),POINTER :: qu(:,:)
162	REAL(rstd),POINTER :: qv(:,:)
163
164	! temporary shared variable
165	REAL(rstd),POINTER :: theta(:,:)
166	REAL(rstd),POINTER :: pk(:,:)
167	REAL(rstd),POINTER :: geopot(:,:)
168	REAL(rstd),POINTER :: convm(:,:)
169	REAL(rstd),POINTER :: wwuu(:,:)
170
171	INTEGER :: ind
172	LOGICAL,SAVE :: first=.TRUE.
173	!$OMP THREADPRIVATE(first)
174
175	IF (first) THEN
176	first=.FALSE.
177	IF(caldyn_eta==eta_mass) THEN
178	CALL init_message(f_ps,req_i1,req_ps)
179	ELSE
180	CALL init_message(f_mass,req_i1,req_mass)
181	END IF
182	CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
183	CALL init_message(f_u,req_e1_vect,req_u)
184	CALL init_message(f_qu,req_e1_scal,req_qu)
185	! Overlapping com/compute (deactivated) : MPI messages need to be sent at first call to caldyn
186	! This is needed only once : the next ones will be sent by timeloop
187	! IF(caldyn_eta==eta_mass) THEN
188	! CALL send_message(f_ps,req_ps)
189	! CALL wait_message(req_ps)
190	! ELSE
191	! CALL send_message(f_mass,req_mass)
192	! CALL wait_message(req_mass)
193	! END IF
194	ENDIF
195
196	CALL trace_start("caldyn")
197
198	IF(caldyn_eta==eta_mass) THEN
199	CALL send_message(f_ps,req_ps) ! COM00
200	ELSE
201	CALL send_message(f_mass,req_mass) ! COM00
202	END IF
203
204	CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
205	CALL send_message(f_u,req_u) ! COM02
206
207	SELECT CASE(caldyn_conserv)
208	CASE(energy) ! energy-conserving
209	DO ind=1,ndomain
210	IF (.NOT. assigned_domain(ind)) CYCLE
211	CALL swap_dimensions(ind)
212	CALL swap_geometry(ind)
213	ps=f_ps(ind)
214	u=f_u(ind)
215	theta_rhodz = f_theta_rhodz(ind)
216	mass=f_mass(ind)
217	theta = f_theta(ind)
218	qu=f_qu(ind)
219	qv=f_qv(ind)
220	CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02
221	ENDDO
222	! CALL checksum(f_mass)
223	! CALL checksum(f_theta)
224
225	CALL send_message(f_qu,req_qu) ! COM03 wait_message in caldyn_horiz
226	! CALL wait_message(req_qu)
227
228	DO ind=1,ndomain
229	IF (.NOT. assigned_domain(ind)) CYCLE
230	CALL swap_dimensions(ind)
231	CALL swap_geometry(ind)
232	ps=f_ps(ind)
233	u=f_u(ind)
234	mass=f_mass(ind)
235	theta = f_theta(ind)
236	qu=f_qu(ind)
237	pk = f_pk(ind)
238	geopot = f_geopot(ind)
239	CALL compute_geopot(ps,mass,theta, pk,geopot)
240	hflux=f_hflux(ind)
241	convm = f_dmass(ind)
242	dtheta_rhodz=f_dtheta_rhodz(ind)
243	du=f_du(ind)
244	CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
245	IF(caldyn_eta==eta_mass) THEN
246	wflux=f_wflux(ind)
247	wwuu=f_wwuu(ind)
248	dps=f_dps(ind)
249	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz(:,:,1), du)
250	END IF
251	ENDDO
252
253	! CALL checksum(f_geopot)
254	! CALL checksum(f_dmass)
255	! CALL checksum(f_pk)
256	! CALL checksum(f_pk)
257
258	CASE(enstrophy) ! enstrophy-conserving
259	DO ind=1,ndomain
260	IF (.NOT. assigned_domain(ind)) CYCLE
261	CALL swap_dimensions(ind)
262	CALL swap_geometry(ind)
263	ps=f_ps(ind)
264	u=f_u(ind)
265	theta_rhodz=f_theta_rhodz(ind)
266	mass=f_mass(ind)
267	theta = f_theta(ind)
268	qu=f_qu(ind)
269	qv=f_qv(ind)
270	CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv)
271	pk = f_pk(ind)
272	geopot = f_geopot(ind)
273	CALL compute_geopot(ps,mass,theta, pk,geopot)
274	hflux=f_hflux(ind)
275	convm = f_dmass(ind)
276	dtheta_rhodz=f_dtheta_rhodz(ind)
277	du=f_du(ind)
278	CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
279	IF(caldyn_eta==eta_mass) THEN
280	wflux=f_wflux(ind)
281	wwuu=f_wwuu(ind)
282	dps=f_dps(ind)
283	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
284	END IF
285	ENDDO
286
287	CASE DEFAULT
288	STOP
289	END SELECT
290
291	!$OMP BARRIER
292	! CALL check_mass_conservation(f_ps,f_dps)
293	CALL trace_end("caldyn")
294	!!$OMP BARRIER
295
296	END SUBROUTINE caldyn
297
298	!-------------------------------- Diagnostics ----------------------------
299
300	SUBROUTINE check_mass_conservation(f_ps,f_dps)
301	USE icosa
302	USE mpipara
303	IMPLICIT NONE
304	TYPE(t_field),POINTER :: f_ps(:)
305	TYPE(t_field),POINTER :: f_dps(:)
306	REAL(rstd),POINTER :: ps(:)
307	REAL(rstd),POINTER :: dps(:)
308	REAL(rstd) :: mass_tot,dmass_tot
309	INTEGER :: ind,i,j,ij
310
311	mass_tot=0
312	dmass_tot=0
313
314	CALL transfert_request(f_dps,req_i1)
315	CALL transfert_request(f_ps,req_i1)
316
317	DO ind=1,ndomain
318	CALL swap_dimensions(ind)
319	CALL swap_geometry(ind)
320
321	ps=f_ps(ind)
322	dps=f_dps(ind)
323
324	DO j=jj_begin,jj_end
325	DO i=ii_begin,ii_end
326	ij=(j-1)*iim+i
327	IF (domain(ind)%own(i,j)) THEN
328	mass_tot=mass_tot+ps(ij)*Ai(ij)/g
329	dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
330	ENDIF
331	ENDDO
332	ENDDO
333
334	ENDDO
335	IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot," dmass_tot ",dmass_tot
336
337	END SUBROUTINE check_mass_conservation
338
339	END MODULE caldyn_gcm_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats: