Context Navigation

source: codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90 @ 899

Last change on this file since 899 was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

Fixed iso c bindings
fixed warnings with -Wall -Wno-aliasing -Wno-unused -Wno-unused-dummy-argument -Wno-maybe-uninitialized -Wno-tabs warnings
Removed all unused variables (-Wunused-variable)
vector%dot_product is now dot_product_3d to avoid compilation warning "dot_product shadows intrinsic" with GCC

File size: 12.3 KB

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

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

Download in other formats: