Context Navigation

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

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

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

Download in other formats: