Context Navigation

interface_icosa_lmdz.f90 @ 4121

Last change on this file since 4121 was 4121, checked in by jgipsl, 6 years ago
Put back bug on mass flux corrected in changeset [310]. This is done only to optain a revision corresponding to what has been used for simulations during the Grand Challenge at Irene/TGCC summer 2018. The bug will be corrected again in next comming commit.
File size: 27.2 KB

Line
1	MODULE interface_icosa_lmdz_mod
2
3	USE field_mod, ONLY: t_field
4	USE transfert_mod, ONLY: t_message
5
6
7	TYPE(t_message),SAVE :: req_u
8	TYPE(t_message),SAVE :: req_dps0, req_dulon0, req_dulat0, req_dTemp0, req_dq0
9
10	TYPE(t_field),POINTER,SAVE :: f_p(:)
11	TYPE(t_field),POINTER,SAVE :: f_pks(:)
12	TYPE(t_field),POINTER,SAVE :: f_pk(:)
13	TYPE(t_field),POINTER,SAVE :: f_p_layer(:)
14	TYPE(t_field),POINTER,SAVE :: f_theta(:)
15	TYPE(t_field),POINTER,SAVE :: f_phi(:)
16	TYPE(t_field),POINTER,SAVE :: f_Temp(:)
17	TYPE(t_field),POINTER,SAVE :: f_ulon(:)
18	TYPE(t_field),POINTER,SAVE :: f_ulat(:)
19	TYPE(t_field),POINTER,SAVE :: f_dulon(:)
20	TYPE(t_field),POINTER,SAVE :: f_dulat(:)
21	TYPE(t_field),POINTER,SAVE :: f_dTemp(:)
22	TYPE(t_field),POINTER,SAVE :: f_dq(:)
23	TYPE(t_field),POINTER,SAVE :: f_dps(:)
24	TYPE(t_field),POINTER,SAVE :: f_duc(:)
25	TYPE(t_field),POINTER,SAVE :: f_bounds_lon(:)
26	TYPE(t_field),POINTER,SAVE :: f_bounds_lat(:)
27
28	INTEGER :: start_clock
29	INTEGER :: stop_clock
30	INTEGER :: count_clock=0
31
32	INTEGER,SAVE :: nbp_phys
33	INTEGER,SAVE :: nbp_phys_glo
34
35
36	CONTAINS
37
38	SUBROUTINE initialize_physics
39	USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
40	! from dynamico
41	USE domain_mod
42	USE dimensions
43	USE mpi_mod
44	USE mpipara
45	USE disvert_mod
46	USE xios_mod
47	USE time_mod , init_time_icosa=> init_time
48	USE transfert_mod
49
50	! from LMDZ
51	USE mod_grid_phy_lmdz, ONLY : unstructured
52	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
53	USE transfert_mod
54	USE physics_distribution_mod, ONLY : init_physics_distribution
55
56
57	IMPLICIT NONE
58	INTEGER :: ind,i,j,ij,pos
59	REAL(rstd),POINTER :: bounds_lon(:,:)
60	REAL(rstd),POINTER :: bounds_lat(:,:)
61
62	REAL(rstd),ALLOCATABLE :: latfi(:)
63	REAL(rstd),ALLOCATABLE :: lonfi(:)
64	REAL(rstd),ALLOCATABLE :: airefi(:)
65	REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
66	REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
67	! REAL(rstd) :: pseudoalt(llm)
68
69	INTEGER :: nbp_phys, nbp_phys_glo
70
71	!$OMP PARALLEL
72	CALL allocate_field(f_bounds_lon,field_t,type_real,6)
73	CALL allocate_field(f_bounds_lat,field_t,type_real,6)
74	CALL allocate_field(f_p,field_t,type_real,llm+1,name="p_in")
75	CALL allocate_field(f_pks,field_t,type_real)
76	CALL allocate_field(f_pk,field_t,type_real,llm)
77	CALL allocate_field(f_p_layer,field_t,type_real,llm,name="p_layer_in")
78	CALL allocate_field(f_theta,field_t,type_real,llm)
79	CALL allocate_field(f_phi,field_t,type_real,llm,name="phi_in")
80	CALL allocate_field(f_Temp,field_t,type_real,llm,name="Temp_in")
81	CALL allocate_field(f_ulon,field_t,type_real,llm,name="ulon_in")
82	CALL allocate_field(f_ulat,field_t,type_real,llm,name="ulat_in")
83	CALL allocate_field(f_dulon,field_t,type_real,llm,name="dulon_out")
84	CALL allocate_field(f_dulat,field_t,type_real,llm,name="dulat_out")
85	CALL allocate_field(f_dTemp,field_t,type_real,llm,name="dTemp_out")
86	CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,name="dq_out")
87	CALL allocate_field(f_dps,field_t,type_real,name="dps_out")
88	CALL allocate_field(f_duc,field_t,type_real,3,llm)
89
90	CALL init_message(f_dps,req_i0,req_dps0)
91	CALL init_message(f_dulon,req_i0,req_dulon0)
92	CALL init_message(f_dulat,req_i0,req_dulat0)
93	CALL init_message(f_dTemp,req_i0,req_dTemp0)
94	CALL init_message(f_dq,req_i0,req_dq0)
95	!$OMP END PARALLEL
96
97	nbp_phys=0
98	DO ind=1,ndomain
99	CALL swap_dimensions(ind)
100	DO j=jj_begin,jj_end
101	DO i=ii_begin,ii_end
102	IF (domain(ind)%own(i,j)) nbp_phys=nbp_phys+1
103	ENDDO
104	ENDDO
105	ENDDO
106
107
108	!initialize LMDZ5 physic mpi decomposition
109	CALL MPI_ALLREDUCE(nbp_phys,nbp_phys_glo,1,MPI_INTEGER,MPI_SUM,comm_icosa,ierr)
110	CALL init_physics_distribution(unstructured, 6, nbp_phys, 1, nbp_phys_glo, llm, comm_icosa)
111
112	DO ind=1,ndomain
113	CALL swap_dimensions(ind)
114	CALL swap_geometry(ind)
115	bounds_lon=f_bounds_lon(ind)
116	bounds_lat=f_bounds_lat(ind)
117	DO j=jj_begin,jj_end
118	DO i=ii_begin,ii_end
119	ij=(j-1)*iim+i
120	CALL xyz2lonlat(xyz_v(ij+z_rup,:), bounds_lon(ij,1), bounds_lat(ij,1))
121	CALL xyz2lonlat(xyz_v(ij+z_up,:), bounds_lon(ij,2), bounds_lat(ij,2))
122	CALL xyz2lonlat(xyz_v(ij+z_lup,:), bounds_lon(ij,3), bounds_lat(ij,3))
123	CALL xyz2lonlat(xyz_v(ij+z_ldown,:), bounds_lon(ij,4), bounds_lat(ij,4))
124	CALL xyz2lonlat(xyz_v(ij+z_down,:), bounds_lon(ij,5), bounds_lat(ij,5))
125	CALL xyz2lonlat(xyz_v(ij+z_rdown,:), bounds_lon(ij,6), bounds_lat(ij,6))
126	ENDDO
127	ENDDO
128	ENDDO
129
130	!$OMP PARALLEL
131	CALL initialize_physics_omp
132	!$OMP END PARALLEL
133
134	CALL xios_set_context
135
136
137
138
139	END SUBROUTINE initialize_physics
140
141
142	SUBROUTINE initialize_physics_omp
143	USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
144	! from dynamico
145	USE domain_mod
146	USE dimensions
147	USE mpi_mod
148	USE mpipara
149	USE disvert_mod
150	USE earth_const, ONLY: scale_height
151	USE xios_mod
152	USE time_mod , init_time_icosa=> init_time
153	USE omp_para
154
155	! from LMDZ
156	USE mod_grid_phy_lmdz, ONLY : unstructured
157	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
158	USE time_phylmdz_mod, ONLY: init_time_lmdz => init_time
159	USE transfert_mod
160	! USE physics_distribution_mod, ONLY : init_physics_distribution
161	USE geometry_mod, ONLY : init_geometry
162	USE vertical_layers_mod, ONLY : init_vertical_layers
163	USE infotrac_phy, ONLY : init_infotrac_phy
164	USE inifis_mod, ONLY : inifis
165	! USE phyaqua_mod, ONLY : iniaqua
166
167
168	IMPLICIT NONE
169
170
171
172	INTEGER :: ind,i,j,k,ij,pos
173	REAL(rstd),POINTER :: bounds_lon(:,:)
174	REAL(rstd),POINTER :: bounds_lat(:,:)
175
176	REAL(rstd),ALLOCATABLE :: latfi(:)
177	REAL(rstd),ALLOCATABLE :: lonfi(:)
178	REAL(rstd),ALLOCATABLE :: airefi(:)
179	REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
180	REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
181	REAL(rstd),ALLOCATABLE :: ind_cell_glo(:)
182
183	REAL(rstd) :: pseudoalt(llm)
184	REAL(rstd) :: aps(llm)
185	REAL(rstd) :: bps(llm)
186	REAL(rstd) :: scaleheight
187
188	INTEGER :: run_length
189	REAL :: day_length ! length of a day (s) ! SAVEd to be OpenMP shared <--- NO!!!!
190	INTEGER :: annee_ref
191	INTEGER :: day_ref
192	INTEGER :: day_ini
193	REAL :: start_time
194	REAL :: physics_timestep
195
196	! Tracer stuff (SAVEd when needed to be OpenMP shared)
197	INTEGER :: nq
198	INTEGER :: nqo, nbtr
199	CHARACTER(len=4) :: type_trac
200	CHARACTER(len=20),ALLOCATABLE :: tname(:) ! tracer short name for restart and diagnostics
201	CHARACTER(len=23),ALLOCATABLE :: ttext(:) ! tracer long name for diagnostics
202	INTEGER,ALLOCATABLE :: niadv(:) ! equivalent dyn / physique
203	INTEGER,ALLOCATABLE :: conv_flg(:) ! conv_flg(it)=0 : convection desactivated for tracer number it
204	INTEGER,ALLOCATABLE :: pbl_flg(:) ! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it
205	CHARACTER(len=8),ALLOCATABLE :: solsym(:) ! tracer name from inca
206	! Isotopes
207	INTEGER,ALLOCATABLE :: nqfils(:)
208	INTEGER,ALLOCATABLE :: nqdesc(:)
209	INTEGER :: nqdesc_tot
210	INTEGER,ALLOCATABLE :: iqfils(:,:)
211	INTEGER,ALLOCATABLE :: iqpere(:)
212	LOGICAL :: ok_isotopes
213	LOGICAL :: ok_iso_verif
214	LOGICAL :: ok_isotrac
215	LOGICAL :: ok_init_iso
216	INTEGER :: niso_possibles
217	REAL,ALLOCATABLE :: tnat(:)
218	REAL,ALLOCATABLE :: alpha_ideal(:)
219	LOGICAL,ALLOCATABLE :: use_iso(:)
220	INTEGER,ALLOCATABLE :: iqiso(:,:)
221	INTEGER,ALLOCATABLE :: iso_num(:)
222	INTEGER,ALLOCATABLE :: iso_indnum(:)
223	INTEGER,ALLOCATABLE :: zone_num(:)
224	INTEGER,ALLOCATABLE :: phase_num(:)
225	INTEGER,ALLOCATABLE :: indnum_fn_num(:)
226	INTEGER,ALLOCATABLE :: index_trac(:,:)
227	INTEGER :: niso
228	INTEGER :: ntraceurs_zone
229	INTEGER :: ntraciso
230
231	TYPE(t_field),POINTER,SAVE :: f_ind_cell_glo(:)
232
233	INTEGER :: iflag_phys
234
235
236	CALL init_distrib_icosa_lmdz
237
238	ALLOCATE(latfi(klon_omp))
239	ALLOCATE(lonfi(klon_omp))
240	ALLOCATE(airefi(klon_omp))
241	ALLOCATE(bounds_latfi(klon_omp,6))
242	ALLOCATE(bounds_lonfi(klon_omp,6))
243	ALLOCATE(ind_cell_glo(klon_omp))
244
245	CALL transfer_icosa_to_lmdz(geom%lat_i,latfi)
246	CALL transfer_icosa_to_lmdz(geom%lon_i,lonfi)
247	CALL transfer_icosa_to_lmdz(f_bounds_lat,bounds_latfi)
248	CALL transfer_icosa_to_lmdz(f_bounds_lon,bounds_lonfi)
249	CALL transfer_icosa_to_lmdz(geom%Ai,airefi)
250
251	CALL allocate_field(f_ind_cell_glo,field_t,type_real)
252
253	DO ind=1,ndomain
254	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
255	CALL swap_dimensions(ind)
256	CALL swap_geometry(ind)
257	DO j=jj_begin,jj_end
258	DO i=ii_begin,ii_end
259	ij=(j-1)*iim+i
260	f_ind_cell_glo(ind)%rval2d(ij)=domain(ind)%assign_cell_glo(i,j)
261	ENDDO
262	ENDDO
263	ENDDO
264
265
266	CALL transfer_icosa_to_lmdz(f_ind_cell_glo,ind_cell_glo)
267	CALL deallocate_field(f_ind_cell_glo)
268
269
270	CALL init_geometry(klon_omp,lonfi, latfi, bounds_lonfi, bounds_latfi, airefi, INT(ind_cell_glo))
271
272	scaleheight=scale_height/1000. ! Atmospheric scale height (km)
273	aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
274	bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
275	pseudoalt(:)=-scaleheight*log(presnivs(:)/preff)
276	CALL init_vertical_layers(llm,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
277
278
279	! Initialize tracer names, numbers, etc. for physics
280	!Config Key = type_trac
281	!Config Desc = Choix de couplage avec model de chimie INCA ou REPROBUS
282	!Config Def = lmdz
283	!Config Help =
284	!Config 'lmdz' = pas de couplage, pur LMDZ
285	!Config 'inca' = model de chime INCA
286	!Config 'repr' = model de chime REPROBUS
287	type_trac = 'lmdz'
288	CALL getin('type_trac',type_trac)
289
290	! allocate some of the tracer arrays
291	ALLOCATE(tname(nqtot))
292	ALLOCATE(ttext(nqtot))
293	ALLOCATE(niadv(nqtot))
294
295	! read "traceur.def" file to know tracer names (and figure out nqo and nbtr)
296	IF (is_mpi_root) THEN
297	!$OMP MASTER
298	OPEN(unit=42,file="traceur.def",form="formatted",status="old",iostat=ierr)
299	IF (ierr==0) THEN
300	READ(42,*) nq ! should be the same as nqtot
301	IF (nq /= nqtot) THEN
302	WRITE(,) "Error: number of tracers in tracer.def should match nqtot!"
303	WRITE(,) " will just use nqtot=",nqtot," tracers"
304	ENDIF
305	DO i=1,nqtot
306	READ(42,*) j,k,tname(i)
307	ENDDO
308	CLOSE(42)
309	ENDIF
310	! figure out how many water tracers are present
311	nqo=0
312	DO i=1,nqtot
313	IF (INDEX(tname(i),"H2O")==1) THEN
314	nqo=nqo+1
315	ENDIF
316	ENDDO
317	nbtr=nqtot-nqo
318	!$OMP END MASTER
319	ENDIF
320
321	CALL bcast(nqo)
322	CALL bcast(nbtr)
323	DO i=1,nqtot
324	CALL bcast(tname(i))
325	ENDDO
326
327	ALLOCATE(conv_flg(nbtr))
328	ALLOCATE(pbl_flg(nbtr))
329	ALLOCATE(solsym(nbtr))
330
331	conv_flg(:) = 1 ! convection activated for all tracers
332	pbl_flg(:) = 1 ! boundary layer activated for all tracers
333	! tracer long names:
334	ttext(1)=trim(tname(1))//"VLH" !'H2OvVLH'
335	DO i=2,nqtot
336	ttext(i)=trim(tname(1))//"VL1"
337	ENDDO
338	solsym(1:nbtr)=tname(nqo+1:nqtot)
339	DO i=1,nqtot
340	niadv(i)=i
341	ENDDO
342	! isotopes
343	ALLOCATE(nqfils(nqtot)) ; nqfils(:)=0
344	ALLOCATE(nqdesc(nqtot)) ; nqdesc(:)=0
345	nqdesc_tot=0
346	ALLOCATE(iqfils(nqtot,nqtot)) ; iqfils(:,:)=0
347	ALLOCATE(iqpere(nqtot)) ; iqpere(:)=0
348	ok_isotopes=.false.
349	ok_iso_verif=.false.
350	ok_isotrac=.false.
351	ok_init_iso=.false.
352	niso_possibles=5
353	niso=0
354	ntraceurs_zone=0
355	ntraciso=0
356	ALLOCATE(tnat(niso_possibles)) ; tnat(:)=0
357	ALLOCATE(alpha_ideal(niso_possibles)) ; alpha_ideal(:)=0
358	ALLOCATE(use_iso(niso_possibles)) ; use_iso(:)=.false.
359	ALLOCATE(iqiso(ntraciso,nqo)) ; iqiso(:,:)=0
360	ALLOCATE(iso_num(nqtot)) ; iso_num(:)=0
361	ALLOCATE(iso_indnum(nqtot)) ; iso_indnum(:)=0
362	ALLOCATE(zone_num(nqtot)) ; zone_num(:)=0
363	ALLOCATE(phase_num(nqtot)) ; phase_num(:)=0
364	ALLOCATE(indnum_fn_num(niso_possibles)) ; indnum_fn_num(:)=0
365	ALLOCATE(index_trac(ntraceurs_zone,niso)) ; index_trac(:,:)=0
366
367	CALL init_infotrac_phy(nqtot,nqo,nbtr,tname,ttext,type_trac,&
368	niadv,conv_flg,pbl_flg,solsym, &
369	nqfils,nqdesc,nqdesc_tot,iqfils,iqpere, &
370	ok_isotopes,ok_iso_verif,ok_isotrac, &
371	ok_init_iso,niso_possibles,tnat, &
372	alpha_ideal,use_iso,iqiso,iso_num, &
373	iso_indnum,zone_num,phase_num, &
374	indnum_fn_num,index_trac, &
375	niso,ntraceurs_zone,ntraciso)
376
377	! Initialize physical constant
378	day_length=86400
379	CALL getin('day_length',day_length)
380	CALL inifis(day_length,radius,g,kappa*cpp,cpp)
381
382
383
384	! init time
385	annee_ref=2015
386	CALL getin("anneeref",annee_ref)
387
388	day_ref=1
389	CALL getin("dayref",day_ref)
390
391	physics_timestep=dt*itau_physics
392	run_length=itaumax*dt
393	ndays=NINT(run_length/day_length)
394
395	day_ini=INT(itau0*dt/day_length)+day_ref
396	start_time= itau0dt/day_length-INT(itau0dt/day_length)
397
398	CALL init_time_lmdz(annee_ref, day_ref, day_ini, start_time, ndays, physics_timestep)
399
400	! Additional initializations for aquaplanets
401	! CALL getin("iflag_phys",iflag_phys)
402	! IF (iflag_phys>=100) THEN
403	! CALL iniaqua(klon_omp, iflag_phys)
404	! END IF
405
406
407	END SUBROUTINE initialize_physics_omp
408
409
410
411
412	SUBROUTINE physics
413	USE icosa
414	USE time_mod
415	USE disvert_mod
416	USE transfert_mod
417	USE mpipara
418	USE xios_mod
419	USE wxios
420	USE trace
421	USE distrib_icosa_lmdz_mod, ONLY : transfer_icosa_to_lmdz, transfer_lmdz_to_icosa
422	USE physics_external_mod, ONLY : it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q
423	USE write_field_mod
424	USE checksum_mod
425	! from LMDZ
426	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
427	USE geometry_mod, ONLY : cell_area
428	USE physiq_mod, ONLY: physiq
429	IMPLICIT NONE
430
431	REAL(rstd),POINTER :: phis(:)
432	REAL(rstd),POINTER :: ps(:)
433	REAL(rstd),POINTER :: theta_rhodz(:,:,:)
434	REAL(rstd),POINTER :: u(:,:)
435	REAL(rstd),POINTER :: wflux(:,:)
436	REAL(rstd),POINTER :: q(:,:,:)
437	REAL(rstd),POINTER :: p(:,:)
438	REAL(rstd),POINTER :: pks(:)
439	REAL(rstd),POINTER :: pk(:,:)
440	REAL(rstd),POINTER :: p_layer(:,:)
441	REAL(rstd),POINTER :: theta(:,:)
442	REAL(rstd),POINTER :: phi(:,:)
443	REAL(rstd),POINTER :: Temp(:,:)
444	REAL(rstd),POINTER :: ulon(:,:)
445	REAL(rstd),POINTER :: ulat(:,:)
446	REAL(rstd),POINTER :: dulon(:,:)
447	REAL(rstd),POINTER :: dulat(:,:)
448	REAL(rstd),POINTER :: dTemp(:,:)
449	REAL(rstd),POINTER :: dq(:,:,:)
450	REAL(rstd),POINTER :: dps(:)
451	REAL(rstd),POINTER :: duc(:,:,:)
452
453
454	INTEGER :: ind,l
455
456	REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
457	!$OMP THREADPRIVATE(ps_phy)
458	REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
459	!$OMP THREADPRIVATE(p_phy)
460	REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
461	!$OMP THREADPRIVATE(p_layer_phy)
462	REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
463	!$OMP THREADPRIVATE(Temp_phy)
464	REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
465	!$OMP THREADPRIVATE(phis_phy)
466	REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
467	!$OMP THREADPRIVATE(phi_phy)
468	REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
469	!$OMP THREADPRIVATE(ulon_phy)
470	REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
471	!$OMP THREADPRIVATE(ulat_phy)
472	REAL(rstd),ALLOCATABLE,SAVE :: rot_phy(:,:)
473	!$OMP THREADPRIVATE(rot_phy)
474	REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
475	!$OMP THREADPRIVATE(q_phy)
476	REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
477	!$OMP THREADPRIVATE(wflux_phy)
478	REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
479	!$OMP THREADPRIVATE(dulon_phy)
480	REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
481	!$OMP THREADPRIVATE(dulat_phy)
482	REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
483	!$OMP THREADPRIVATE(dTemp_phy)
484	REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
485	!$OMP THREADPRIVATE(dq_phy)
486	REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
487	!$OMP THREADPRIVATE(dps_phy)
488	REAL(rstd) :: dtphy
489	LOGICAL :: debut
490	LOGICAL :: lafin
491	LOGICAL,SAVE :: first=.TRUE.
492	!$OMP THREADPRIVATE(first)
493
494
495	IF(first) THEN
496	debut=.TRUE.
497	ELSE
498	debut=.FALSE.
499	ENDIF
500
501
502	IF(it-itau0>=itaumax) THEN
503	lafin=.TRUE.
504	ELSE
505	lafin=.FALSE.
506	ENDIF
507
508	IF (first) THEN
509	first=.FALSE.
510	CALL init_message(f_u,req_e1_vect,req_u)
511	ALLOCATE(ps_phy(klon_omp))
512	ALLOCATE(p_phy(klon_omp,llm+1))
513	ALLOCATE(p_layer_phy(klon_omp,llm))
514	ALLOCATE(Temp_phy(klon_omp,llm))
515	ALLOCATE(phis_phy(klon_omp))
516	ALLOCATE(phi_phy(klon_omp,llm))
517	ALLOCATE(ulon_phy(klon_omp,llm))
518	ALLOCATE(ulat_phy(klon_omp,llm))
519	ALLOCATE(rot_phy(klon_omp,llm))
520	ALLOCATE(q_phy(klon_omp,llm,nqtot))
521	ALLOCATE(wflux_phy(klon_omp,llm))
522	ALLOCATE(dulon_phy(klon_omp,llm))
523	ALLOCATE(dulat_phy(klon_omp,llm))
524	ALLOCATE(dTemp_phy(klon_omp,llm))
525	ALLOCATE(dq_phy(klon_omp,llm,nqtot))
526	ALLOCATE(dps_phy(klon_omp))
527	!$OMP BARRIER
528	ENDIF
529
530
531	!$OMP MASTER
532	! CALL update_calendar(it)
533	!$OMP END MASTER
534	!$OMP BARRIER
535	dtphy=itau_physics*dt
536
537
538
539	CALL transfert_message(f_u,req_u)
540
541	DO ind=1,ndomain
542	CALL swap_dimensions(ind)
543	IF (assigned_domain(ind)) THEN
544	CALL swap_geometry(ind)
545
546	phis=f_phis(ind)
547	ps=f_ps(ind)
548	theta_rhodz=f_theta_rhodz(ind)
549	u=f_u(ind)
550	q=f_q(ind)
551	wflux=f_wflux(ind)
552	p=f_p(ind)
553	pks=f_pks(ind)
554	pk=f_pk(ind)
555	p_layer=f_p_layer(ind)
556	theta=f_theta(ind)
557	phi=f_phi(ind)
558	Temp=f_Temp(ind)
559	ulon=f_ulon(ind)
560	ulat=f_ulat(ind)
561
562	CALL grid_icosa_to_physics
563
564	ENDIF
565	ENDDO
566
567	!$OMP BARRIER
568	!$OMP MASTER
569	CALL SYSTEM_CLOCK(start_clock)
570	!$OMP END MASTER
571	CALL trace_start("physic")
572	! CALL trace_off()
573
574
575	! CALL writeField("p_in",f_p)
576	! CALL writeField("p_layer_in",f_p_layer)
577	! CALL writeField("phi_in",f_phi)
578	! CALL writeField("phis_in",f_phis)
579	! CALL writeField("ulon_in",f_ulon)
580	! CALL writeField("ulat_in",f_ulat)
581	! CALL writeField("Temp_in",f_Temp)
582	! CALL writeField("q_in",f_q)
583	! CALL writeField("wflux_in",f_wflux)
584
585	! CALL checksum(f_p)
586	! CALL checksum(f_p_layer)
587	! CALL checksum(f_phi)
588	! CALL checksum(f_phis)
589	! CALL checksum(f_ulon)
590	! CALL checksum(f_ulat)
591	! CALL checksum(f_Temp)
592	! CALL checksum(f_q)
593	! CALL checksum(f_wflux)
594
595	CALL transfer_icosa_to_lmdz(f_p , p_phy)
596	CALL transfer_icosa_to_lmdz(f_p_layer, p_layer_phy)
597	CALL transfer_icosa_to_lmdz(f_phi , phi_phy)
598	CALL transfer_icosa_to_lmdz(f_phis , phis_phy )
599	CALL transfer_icosa_to_lmdz(f_ulon , ulon_phy )
600	CALL transfer_icosa_to_lmdz(f_ulat , ulat_phy)
601	CALL transfer_icosa_to_lmdz(f_Temp , Temp_phy)
602	CALL transfer_icosa_to_lmdz(f_q , q_phy)
603	CALL transfer_icosa_to_lmdz(f_wflux , wflux_phy)
604
605	DO l=1,llm
606	wflux_phy(:,l)=wflux_phy(:,l)*cell_area(:)
607	phi_phy(:,l)=phi_phy(:,l)-phis_phy(:)
608	ENDDO
609
610	CALL wxios_set_context()
611
612	! Ehouarn: rot_phy() not implemented!! Set it to zero for now
613	rot_phy(:,:)=0
614	CALL physiq(klon_omp, llm, debut, lafin, dtphy, &
615	p_phy, p_layer_phy, phi_phy, phis_phy, presnivs, &
616	ulon_phy, ulat_phy, rot_phy, Temp_phy, q_phy, wflux_phy, &
617	dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy)
618
619	CALL transfer_lmdz_to_icosa(dulon_phy, f_dulon )
620	CALL transfer_lmdz_to_icosa(dulat_phy, f_dulat )
621	CALL transfer_lmdz_to_icosa(dTemp_phy, f_dTemp )
622	CALL transfer_lmdz_to_icosa(dq_phy , f_dq )
623	CALL transfer_lmdz_to_icosa(dps_phy , f_dps )
624
625	! CALL writeField("dulon_out",f_dulon)
626	! CALL writeField("dulat_out",f_dulat)
627	! CALL writeField("dTemp_out",f_dTemp)
628	! CALL writeField("dq_out",f_dq)
629	! CALL writeField("dps_out",f_dps)
630
631	! CALL checksum(f_dulon)
632	! CALL checksum(f_dulat)
633	! CALL checksum(f_dTemp)
634	! CALL checksum(f_dq)
635	! CALL checksum(f_dps)
636
637	CALL send_message(f_dps,req_dps0)
638	CALL send_message(f_dulon,req_dulon0)
639	CALL send_message(f_dulat,req_dulat0)
640	CALL send_message(f_dTemp,req_dTemp0)
641	CALL send_message(f_dq,req_dq0)
642
643	CALL wait_message(req_dps0)
644	CALL wait_message(req_dulon0)
645	CALL wait_message(req_dulat0)
646	CALL wait_message(req_dTemp0)
647	CALL wait_message(req_dq0)
648
649
650	! CALL trace_on()
651	CALL trace_end("physic")
652	!$OMP MASTER
653	CALL SYSTEM_CLOCK(stop_clock)
654	count_clock=count_clock+stop_clock-start_clock
655	!$OMP END MASTER
656
657	!$OMP BARRIER
658
659	DO ind=1,ndomain
660	CALL swap_dimensions(ind)
661	IF (assigned_domain(ind)) THEN
662	CALL swap_geometry(ind)
663
664	theta_rhodz=f_theta_rhodz(ind)
665	u=f_u(ind)
666	q=f_q(ind)
667	ps=f_ps(ind)
668	dulon=f_dulon(ind)
669	dulat=f_dulat(ind)
670	Temp=f_temp(ind)
671	dTemp=f_dTemp(ind)
672	dq=f_dq(ind)
673	dps=f_dps(ind)
674	duc=f_duc(ind)
675	p=f_p(ind)
676	pks=f_pks(ind)
677	pk=f_pk(ind)
678
679	CALL grid_physics_to_icosa
680	ENDIF
681	ENDDO
682
683	!$OMP BARRIER
684	CALL xios_set_context
685
686
687	CONTAINS
688
689	SUBROUTINE grid_icosa_to_physics
690	USE pression_mod
691	USE exner_mod
692	USE theta2theta_rhodz_mod
693	USE geopotential_mod
694	USE wind_mod
695	USE omp_para
696	IMPLICIT NONE
697
698	REAL(rstd) :: uc(3)
699	INTEGER :: i,j,ij,l
700
701
702	! compute pression
703
704	DO l = ll_begin,ll_endp1
705	DO j=jj_begin,jj_end
706	DO i=ii_begin,ii_end
707	ij=(j-1)*iim+i
708	p(ij,l) = ap(l) + bp(l) * ps(ij)
709	ENDDO
710	ENDDO
711	ENDDO
712
713	!$OMP BARRIER
714
715	! compute exner
716
717	IF (is_omp_first_level) THEN
718	DO j=jj_begin,jj_end
719	DO i=ii_begin,ii_end
720	ij=(j-1)*iim+i
721	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
722	ENDDO
723	ENDDO
724	ENDIF
725
726	! 3D : pk
727	DO l = ll_begin,ll_end
728	DO j=jj_begin,jj_end
729	DO i=ii_begin,ii_end
730	ij=(j-1)*iim+i
731	pk(ij,l) = cpp * ((.5/preff)(p(ij,l)+p(ij,l+1))) * kappa
732	ENDDO
733	ENDDO
734	ENDDO
735
736	!$OMP BARRIER
737
738	! compute theta, temperature and pression at layer
739	DO l = ll_begin, ll_end
740	DO j=jj_begin,jj_end
741	DO i=ii_begin,ii_end
742	ij=(j-1)*iim+i
743	theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g)
744	Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp
745	p_layer(ij,l)=preff(pk(ij,l)/cpp)*(1./kappa)
746	ENDDO
747	ENDDO
748	ENDDO
749
750
751	!!! Compute geopotential
752
753	! for first layer
754	IF (is_omp_first_level) THEN
755	DO j=jj_begin,jj_end
756	DO i=ii_begin,ii_end
757	ij=(j-1)*iim+i
758	phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
759	ENDDO
760	ENDDO
761	ENDIF
762	!!-> implicit flush on phi(:,1)
763
764	! for other layers
765	DO l = ll_beginp1, ll_end
766	DO j=jj_begin,jj_end
767	DO i=ii_begin,ii_end
768	ij=(j-1)*iim+i
769	phi(ij,l) = 0.5 * ( theta(ij,l) + theta(ij,l-1) ) &
770	* ( pk(ij,l-1) - pk(ij,l) )
771	ENDDO
772	ENDDO
773	ENDDO
774
775	!$OMP BARRIER
776
777
778	IF (is_omp_first_level) THEN
779	DO l = 2, llm
780	DO j=jj_begin,jj_end
781	! ---> Bug compilo intel ici en openmp
782	! ---> Couper la boucle
783	IF (j==jj_end+1) PRINT*,"this message must not be printed"
784	DO i=ii_begin,ii_end
785	ij=(j-1)*iim+i
786	phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
787	ENDDO
788	ENDDO
789	ENDDO
790	! --> IMPLICIT FLUSH on phi --> non
791	ENDIF
792
793	! compute wind centered lon lat compound
794	DO l=ll_begin,ll_end
795	DO j=jj_begin,jj_end
796	DO i=ii_begin,ii_end
797	ij=(j-1)*iim+i
798	uc(:)=1/Ai(ij)* &
799	( ne(ij,right)u(ij+u_right,l)le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:)) &
800	+ ne(ij,rup)u(ij+u_rup,l)le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:)) &
801	+ ne(ij,lup)u(ij+u_lup,l)le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:)) &
802	+ ne(ij,left)u(ij+u_left,l)le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:)) &
803	+ ne(ij,ldown)u(ij+u_ldown,l)le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))&
804	+ ne(ij,rdown)u(ij+u_rdown,l)le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:)))
805	ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
806	ulat(ij,l)=sum(uc(:)*elat_i(ij,:))
807	ENDDO
808	ENDDO
809	ENDDO
810
811	!$OMP BARRIER
812	END SUBROUTINE grid_icosa_to_physics
813
814
815	SUBROUTINE grid_physics_to_icosa
816	USE theta2theta_rhodz_mod
817	USE omp_para
818	IMPLICIT NONE
819	INTEGER :: i,j,ij,l,iq
820
821	DO l=ll_begin,ll_end
822	DO j=jj_begin,jj_end
823	DO i=ii_begin,ii_end
824	ij=(j-1)*iim+i
825	duc(ij,:,l)=dulon(ij,l)elon_i(ij,:)+dulat(ij,l)elat_i(ij,:)
826	ENDDO
827	ENDDO
828	ENDDO
829
830	DO l=ll_begin,ll_end
831	DO j=jj_begin,jj_end
832	DO i=ii_begin,ii_end
833	ij=(j-1)*iim+i
834	u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5(duc(ij,:,l) + duc(ij+t_right,:,l))ep_e(ij+u_right,:) )
835	u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5(duc(ij,:,l) + duc(ij+t_lup,:,l))ep_e(ij+u_lup,:) )
836	u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphysum( 0.5(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
837	ENDDO
838	ENDDO
839	ENDDO
840
841	DO l=ll_begin,ll_end
842	DO j=jj_begin,jj_end
843	DO i=ii_begin,ii_end
844	ij=(j-1)*iim+i
845	Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
846	ENDDO
847	ENDDO
848	ENDDO
849
850	DO iq=1,nqtot
851	DO l=ll_begin,ll_end
852	DO j=jj_begin,jj_end
853	DO i=ii_begin,ii_end
854	ij=(j-1)*iim+i
855	q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
856	ENDDO
857	ENDDO
858	ENDDO
859	ENDDO
860
861	!$OMP BARRIER
862
863	IF (is_omp_first_level) THEN
864	DO j=jj_begin,jj_end
865	DO i=ii_begin,ii_end
866	ij=(j-1)*iim+i
867	ps(ij)=ps(ij)+ dtphy * dps(ij)
868	ENDDO
869	ENDDO
870	ENDIF
871
872	! CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
873
874	! compute pression
875	!$OMP BARRIER
876	DO l = ll_begin,ll_endp1
877	DO j=jj_begin,jj_end
878	DO i=ii_begin,ii_end
879	ij=(j-1)*iim+i
880	p(ij,l) = ap(l) + bp(l) * ps(ij)
881	ENDDO
882	ENDDO
883	ENDDO
884
885	!$OMP BARRIER
886
887	! compute exner
888
889	IF (is_omp_first_level) THEN
890	DO j=jj_begin,jj_end
891	DO i=ii_begin,ii_end
892	ij=(j-1)*iim+i
893	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
894	ENDDO
895	ENDDO
896	ENDIF
897
898	! 3D : pk
899	DO l = ll_begin,ll_end
900	DO j=jj_begin,jj_end
901	DO i=ii_begin,ii_end
902	ij=(j-1)*iim+i
903	pk(ij,l) = cpp * ((.5/preff)(p(ij,l)+p(ij,l+1))) * kappa
904	ENDDO
905	ENDDO
906	ENDDO
907
908	!$OMP BARRIER
909
910	! compute theta, temperature and pression at layer
911	DO l = ll_begin, ll_end
912	DO j=jj_begin,jj_end
913	DO i=ii_begin,ii_end
914	ij=(j-1)*iim+i
915	theta_rhodz(ij,l,1) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
916	ENDDO
917	ENDDO
918	ENDDO
919
920	END SUBROUTINE grid_physics_to_icosa
921
922
923
924	END SUBROUTINE physics
925
926
927
928
929
930	END MODULE interface_icosa_lmdz_mod

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

Download in other formats:

Original Format