Context Navigation

caldyn_gcm.f90 @ 519

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

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

Context Navigation

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

Download in other formats: