Context Navigation

source: CONFIG/UNIFORM/v7/ICOLMDZOR_v7/SOURCES/ICOSA_LMDZ/src/phylmd/interface_icosa_lmdz.f90 @ 5878

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

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

Download in other formats: