Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

interface_icosa_lmdz.F90 @ 5875

Last change on this file since 5875 was 5590, checked in by yann meurdesoif, 8 months ago
Add phylmdiso directory for testing ICOLMDZISO configuration. YM
Property svn:eol-style set to `native`
File size: 39.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	REAL,ALLOCATABLE,SAVE :: q_ave(:, :)
38	!$OMP THREADPRIVATE(q_ave)
39
40	REAL, PARAMETER :: q_epsilon = 1e-8
41
42	CONTAINS
43
44	SUBROUTINE pre_initialize_physics
45	USE etat0_plugin_mod
46	USE icolmdz_etat0, ONLY: init_etat0_lmdz => init_etat0, etat0_lmdz => etat0
47	USE icolmdz_param_gravity_wave, ONLY: init_param_gravity_wave
48	USE isotopes_mod, ONLY : using_iso
49	IMPLICIT NONE
50	!$OMP PARALLEL
51
52	init_etat0_plugin => init_etat0_lmdz
53	etat0_plugin => etat0_lmdz
54
55	CALL init_param_gravity_wave
56	!$OMP END PARALLEL
57	END SUBROUTINE pre_initialize_physics
58
59	SUBROUTINE initialize_physics
60	USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
61	! from dynamico
62	USE domain_mod
63	USE dimensions
64	USE mpi_mod
65	USE mpipara
66	USE disvert_mod
67	USE xios_mod
68	USE time_mod , init_time_icosa=> init_time
69	USE transfert_mod
70	USE nudging_mod, ONLY : lam_halo_scheme
71
72	! from LMDZ
73	USE mod_grid_phy_lmdz, ONLY : unstructured
74	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
75	USE transfert_mod
76	USE physics_distribution_mod, ONLY : init_physics_distribution
77	USE infotrac_phy, ONLY: init_infotrac_phy
78	USE icolmdz_param_gravity_wave, ONLY: param_gravity_wave
79
80
81
82	IMPLICIT NONE
83	INTEGER :: ind,i,j,ij,pos,h
84	REAL(rstd),POINTER :: bounds_lon(:,:)
85	REAL(rstd),POINTER :: bounds_lat(:,:)
86
87	REAL(rstd),ALLOCATABLE :: latfi(:)
88	REAL(rstd),ALLOCATABLE :: lonfi(:)
89	REAL(rstd),ALLOCATABLE :: airefi(:)
90	REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
91	REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
92	LOGICAL ,ALLOCATABLE :: outside(:,:)
93	LOGICAL ,ALLOCATABLE :: outside_tmp(:,:)
94	LOGICAL ,POINTER :: out(:,:)
95	! REAL(rstd) :: pseudoalt(llm)
96
97	INTEGER :: nbp_phys, nbp_phys_glo
98
99	!$OMP PARALLEL
100	CALL allocate_field(f_bounds_lon,field_t,type_real,6)
101	CALL allocate_field(f_bounds_lat,field_t,type_real,6)
102	CALL allocate_field(f_p,field_t,type_real,llm+1,name="p_in")
103	CALL allocate_field(f_pks,field_t,type_real)
104	CALL allocate_field(f_pk,field_t,type_real,llm)
105	CALL allocate_field(f_p_layer,field_t,type_real,llm,name="p_layer_in")
106	CALL allocate_field(f_theta,field_t,type_real,llm)
107	CALL allocate_field(f_phi,field_t,type_real,llm,name="phi_in")
108	CALL allocate_field(f_Temp,field_t,type_real,llm,name="Temp_in")
109	CALL allocate_field(f_ulon,field_t,type_real,llm,name="ulon_in")
110	CALL allocate_field(f_ulat,field_t,type_real,llm,name="ulat_in")
111	CALL allocate_field(f_vort,field_z,type_real,llm,name="vort_in")
112	CALL allocate_field(f_vortc,field_t,type_real,llm,name="vortc_in")
113	CALL allocate_field(f_dulon,field_t,type_real,llm,name="dulon_out")
114	CALL allocate_field(f_dulat,field_t,type_real,llm,name="dulat_out")
115	CALL allocate_field(f_dTemp,field_t,type_real,llm,name="dTemp_out")
116	CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,name="dq_out")
117	CALL allocate_field(f_dps,field_t,type_real,name="dps_out")
118	CALL allocate_field(f_duc,field_t,type_real,3,llm)
119
120	CALL init_message(f_dps,req_i0,req_dps0)
121	CALL init_message(f_dulon,req_i0,req_dulon0)
122	CALL init_message(f_dulat,req_i0,req_dulat0)
123	CALL init_message(f_dTemp,req_i0,req_dTemp0)
124	CALL init_message(f_dq,req_i0,req_dq0)
125	!$OMP END PARALLEL
126
127	nbp_phys=0
128	DO ind=1,ndomain
129	CALL swap_dimensions(ind)
130
131	ALLOCATE(outside(ii_begin:ii_end,jj_begin:jj_end)) ! for limited area : don't take cells arround the border
132	ALLOCATE(outside_tmp(ii_begin-1:ii_end+1,jj_begin-1:jj_end+1)) ! for limited area : don't take cells arround the border
133	out=>domain(ind)%outside
134	DO j=jj_begin,jj_end
135	DO i=ii_begin,ii_end
136	outside(i,j)= out(i+1,j) .OR. & ! right
137	out(i,j+1 ) .OR. & ! rup
138	out(i-1 ,j+1) .OR. & ! lup
139	out(i-1 ,j) .OR. & !left
140	out(i ,j-1) .OR. & !ldown
141	out(i+1,j-1) !rdown
142	ENDDO
143	ENDDO
144
145	outside_tmp=.FALSE.
146	outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside
147
148	DO h=1,lam_halo_scheme-1 ! do not compute physic on limited area halo
149	DO j=jj_begin,jj_end
150	DO i=ii_begin,ii_end
151	outside(i,j) = outside_tmp(i,j) .OR. &
152	outside_tmp(i+1,j) .OR. & ! right
153	outside_tmp(i,j+1 ) .OR. & ! rup
154	outside_tmp(i-1 ,j+1) .OR. & ! lup
155	outside_tmp(i-1 ,j) .OR. & !left
156	outside_tmp(i ,j-1) .OR. & !ldown
157	outside_tmp(i+1,j-1) !rdown
158	ENDDO
159	ENDDO
160	outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside
161	ENDDO
162
163	DO j=jj_begin,jj_end
164	DO i=ii_begin,ii_end
165	IF (domain(ind)%own(i,j) .AND. .NOT.outside(i,j)) nbp_phys=nbp_phys+1
166	ENDDO
167	ENDDO
168	DEALLOCATE(outside)
169	DEALLOCATE(outside_tmp)
170	ENDDO
171
172
173	!initialize LMDZ5 physic mpi decomposition
174	CALL MPI_ALLREDUCE(nbp_phys,nbp_phys_glo,1,MPI_INTEGER,MPI_SUM,comm_icosa,ierr)
175	CALL init_physics_distribution(unstructured, 6, nbp_phys, 1, nbp_phys_glo, llm, comm_icosa)
176
177	DO ind=1,ndomain
178	CALL swap_dimensions(ind)
179	CALL swap_geometry(ind)
180	bounds_lon=f_bounds_lon(ind)
181	bounds_lat=f_bounds_lat(ind)
182	DO j=jj_begin,jj_end
183	DO i=ii_begin,ii_end
184	ij=(j-1)*iim+i
185	CALL xyz2lonlat(xyz_v(ij+z_rup,:), bounds_lon(ij,1), bounds_lat(ij,1))
186	CALL xyz2lonlat(xyz_v(ij+z_up,:), bounds_lon(ij,2), bounds_lat(ij,2))
187	CALL xyz2lonlat(xyz_v(ij+z_lup,:), bounds_lon(ij,3), bounds_lat(ij,3))
188	CALL xyz2lonlat(xyz_v(ij+z_ldown,:), bounds_lon(ij,4), bounds_lat(ij,4))
189	CALL xyz2lonlat(xyz_v(ij+z_down,:), bounds_lon(ij,5), bounds_lat(ij,5))
190	CALL xyz2lonlat(xyz_v(ij+z_rdown,:), bounds_lon(ij,6), bounds_lat(ij,6))
191	ENDDO
192	ENDDO
193	ENDDO
194
195	CALL init_infotrac_phy
196
197
198	!$OMP PARALLEL
199	CALL initialize_physics_omp
200	CALL param_gravity_wave
201	!$OMP END PARALLEL
202
203	CALL xios_set_context
204
205	END SUBROUTINE initialize_physics
206
207
208	SUBROUTINE initialize_physics_omp
209	USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
210	! from dynamico
211	USE domain_mod
212	USE dimensions
213	USE mpi_mod
214	USE mpipara
215	USE disvert_mod
216	USE earth_const, ONLY: scale_height
217	USE xios_mod
218	USE time_mod , init_time_icosa=> init_time
219	USE omp_para
220
221	! from LMDZ
222	USE mod_grid_phy_lmdz, ONLY : unstructured, klon_glo
223	USE mod_phys_lmdz_para, ONLY: klon_omp, reduce_min_lmdz => reduce_min , gather_lmdz => gather , bcast_lmdz => bcast
224	USE time_phylmdz_mod, ONLY: init_time_lmdz => init_time
225	USE transfert_mod
226	! USE physics_distribution_mod, ONLY : init_physics_distribution
227	USE geometry_mod, ONLY : init_geometry
228	USE vertical_layers_mod, ONLY : init_vertical_layers
229	USE infotrac_phy, ONLY : init_infotrac_phy
230	USE inifis_mod, ONLY : inifis
231	USE readTracFiles_mod, ONLY: trac_type, isot_type
232	USE tracer_icosa_mod, ONLY : tracs
233	USE readTracFiles_mod, ONLY: delPhase
234	! USE phyaqua_mod, ONLY : iniaqua
235
236
237	IMPLICIT NONE
238
239
240
241	INTEGER :: ind,i,j,k,ij,pos
242	REAL(rstd),POINTER :: bounds_lon(:,:)
243	REAL(rstd),POINTER :: bounds_lat(:,:)
244
245	REAL(rstd),ALLOCATABLE :: latfi(:)
246	REAL(rstd),ALLOCATABLE :: lonfi(:)
247	REAL(rstd),ALLOCATABLE :: airefi(:)
248	REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
249	REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
250	REAL(rstd),ALLOCATABLE :: ind_cell_glo_r(:)
251	INTEGER, ALLOCATABLE :: ind_cell_glo(:)
252	INTEGER, ALLOCATABLE :: ind_cell_glo_tot(:)
253	INTEGER, ALLOCATABLE :: cell_glo_tot(:)
254	INTEGER :: ncell_glo_tot
255
256	REAL(rstd) :: pseudoalt(llm)
257	REAL(rstd) :: aps(llm)
258	REAL(rstd) :: bps(llm)
259	REAL(rstd) :: scaleheight
260
261	INTEGER :: run_length
262	REAL :: day_length ! length of a day (s) ! SAVEd to be OpenMP shared <--- NO!!!!
263	INTEGER :: annee_ref
264	INTEGER :: day_ref
265	INTEGER :: day_ini
266	REAL :: start_time
267	REAL :: physics_timestep
268
269	! Tracer stuff (SAVEd when needed to be OpenMP shared)
270	INTEGER :: nq
271	INTEGER :: nqo, nbtr, nbtr_inca
272	CHARACTER(len=256) :: type_trac
273	INTEGER,ALLOCATABLE :: conv_flg(:) ! conv_flg(it)=0 : convection desactivated for tracer number it
274	INTEGER,ALLOCATABLE :: pbl_flg(:) ! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it
275	CHARACTER(len=8),ALLOCATABLE :: solsym(:) ! tracer name from inca
276
277	TYPE(t_field),POINTER,SAVE :: f_ind_cell_glo(:)
278
279	INTEGER :: iflag_phys
280
281	INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca ! index of horizontal trasport schema
282	INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca ! index of vertical trasport schema
283
284	TYPE(trac_type) :: tracers_ico2lmdz(nqtot) !=== TRACERS DESCRIPTORS VECTOR
285	TYPE(isot_type) :: isotopes_ico2lmdz(1) !=== ISOTOPES PARAMETERS VECTOR
286	INTEGER :: iq
287
288	CHARACTER(LEN=3) :: descrq(30) !--- Advection scheme description tags
289	logical, save :: first = .TRUE.
290
291
292	CALL init_distrib_icosa_lmdz
293
294	ALLOCATE(latfi(klon_omp))
295	ALLOCATE(lonfi(klon_omp))
296	ALLOCATE(airefi(klon_omp))
297	ALLOCATE(bounds_latfi(klon_omp,6))
298	ALLOCATE(bounds_lonfi(klon_omp,6))
299	ALLOCATE(ind_cell_glo_r(klon_omp))
300	ALLOCATE(ind_cell_glo(klon_omp))
301	ALLOCATE(q_ave(llm,nqtot))
302
303	CALL transfer_icosa_to_lmdz(geom%lat_i,latfi)
304	CALL transfer_icosa_to_lmdz(geom%lon_i,lonfi)
305	CALL transfer_icosa_to_lmdz(f_bounds_lat,bounds_latfi)
306	CALL transfer_icosa_to_lmdz(f_bounds_lon,bounds_lonfi)
307	CALL transfer_icosa_to_lmdz(geom%Ai,airefi)
308
309	CALL allocate_field(f_ind_cell_glo,field_t,type_real)
310
311	DO ind=1,ndomain
312	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
313	CALL swap_dimensions(ind)
314	CALL swap_geometry(ind)
315	DO j=jj_begin,jj_end
316	DO i=ii_begin,ii_end
317	ij=(j-1)*iim+i
318	f_ind_cell_glo(ind)%rval2d(ij)=domain(ind)%assign_cell_glo(i,j)
319	ENDDO
320	ENDDO
321	ENDDO
322
323
324	CALL transfer_icosa_to_lmdz(f_ind_cell_glo,ind_cell_glo_r)
325	CALL deallocate_field(f_ind_cell_glo)
326	ind_cell_glo=INT(ind_cell_glo_r)
327	DEALLOCATE(ind_cell_glo_r)
328
329
330	CALL reduce_min_lmdz(MINVAL(-ind_cell_glo),ncell_glo_tot) ! reduce_max does not exist in lmdz, use reduce_min
331	CALL bcast_lmdz(ncell_glo_tot)
332	ncell_glo_tot=-ncell_glo_tot
333	ALLOCATE(cell_glo_tot(0:ncell_glo_tot))
334	ALLOCATE(ind_cell_glo_tot(klon_glo))
335	CALL gather_lmdz(ind_cell_glo,ind_cell_glo_tot)
336	CALL bcast_lmdz(ind_cell_glo_tot)
337
338	cell_glo_tot=-1
339	DO i=1,klon_glo
340	cell_glo_tot(ind_cell_glo_tot(i))= 0
341	ENDDO
342
343	pos=0
344	DO i=0,ncell_glo_tot
345	IF (cell_glo_tot(i)/=-1) THEN
346	cell_glo_tot(i) = pos
347	pos=pos + 1
348	ENDIF
349	ENDDO
350
351	DO i=1,klon_omp
352	ind_cell_glo(i)=cell_glo_tot(ind_cell_glo(i))
353	ENDDO
354
355	ind_cell_glo = ind_cell_glo + 1 ! lmdz expect global indices begining to 1 not 0
356
357
358	! CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_mpi,1,MPI_INTEGER,comm_icosa,ierr)
359	!
360	! displ(0)=0
361	! DO i=1,mpi_size-1
362	! displ(i)=displ(i-1)+ncell_mpi(i-1)
363	! ENDDO
364	!
365	! ALLOCATE(ind_glo_tot(ncell_tot))
366	! ALLOCATE(cell_glo_tot(0:ncell_glo_tot-1))
367	!
368	! cell_glo_tot(:)= -1
369	! CALL MPI_ALLGATHERV(ind_glo, ncell, MPI_INTEGER, ind_glo_tot, ncell_mpi, displ, MPI_INTEGER, comm_icosa,ierr)
370	!
371	! DO i=1,ncell_tot
372	! cell_glo_tot(ind_glo_tot(i))= 0
373	! ENDDO
374	!
375	! ncell_glo=0
376	! DO i=0,ncell_glo_tot-1
377	! IF (cell_glo_tot(i)/=-1) THEN
378	! cell_glo_tot(i) = ncell_glo
379	! ncell_glo=ncell_glo + 1
380	! ENDIF
381	! ENDDO
382	!
383	! DO i=1,ncell
384	! ind_glo(i)=cell_glo_tot(ind_glo(i))
385	! ENDDO
386
387
388	CALL init_geometry(klon_omp,lonfi, latfi, bounds_lonfi, bounds_latfi, airefi, ind_cell_glo)
389
390	scaleheight=scale_height/1000. ! Atmospheric scale height (km)
391	aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
392	bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
393	pseudoalt(:)=-scaleheight*log(presnivs(:)/preff)
394	CALL init_vertical_layers(llm,preff,scaleheight,ap,bp,aps,bps,presnivs,presinter,pseudoalt)
395
396	! Initialize tracer names, numbers, etc. for physics
397	!Config Key = type_trac
398	!Config Desc = Choix de couplage avec model de chimie INCA ou REPROBUS
399	!Config Def = lmdz
400	!Config Help =
401	!Config 'lmdz' = pas de couplage, pur LMDZ
402	!Config 'lmdz\|inca' = model de chime INCA
403	!Config 'lmdz\|repr' = model de chime REPROBUS
404	type_trac = 'lmdz'
405	CALL getin('type_trac',type_trac)
406
407	descrq( 1: 2) = ['LMV','BAK']
408	descrq(10:20) = ['VL1','VLP','FH1','FH2','VLH',' ','PPM','PPS','PPP',' ','SLP']
409	descrq(30) = 'PRA'
410
411	nqo = 0
412	DO iq=1,nqtot
413
414	tracers_ico2lmdz(iq)%name = tracs(iq)%name
415
416	tracers_ico2lmdz(iq)%gen0Name = tracs(iq)%name
417	tracers_ico2lmdz(iq)%phase = tracs(iq)%phase
418	tracers_ico2lmdz(iq)%iadv = tracs(iq)%iadv
419
420	IF (tracs(iq)%component .eq. "dynamico") then
421	tracers_ico2lmdz(iq)%component='lmdz'
422	ELSE
423	tracers_ico2lmdz(iq)%component=tracs(iq)%component
424	ENDIF
425
426	if (tracers_ico2lmdz(iq)%iadv .ne. 0 ) tracers_ico2lmdz(iq)%isAdvected=.true.
427
428	tracers_ico2lmdz(iq)%longName = tracers_ico2lmdz(iq)%name
429	IF(tracers_ico2lmdz(iq)%iadv > 0) tracers_ico2lmdz(iq)%longName=TRIM(tracers_ico2lmdz(iq)%name)//descrq(tracers_ico2lmdz(iq)%iadv)
430
431	tracers_ico2lmdz(iq)%isInPhysics= delPhase(tracers_ico2lmdz(iq)%gen0Name) /= 'H2O' .OR. tracers_ico2lmdz(iq)%component /= 'lmdz'
432	tracers_ico2lmdz(iq)%iGeneration = 0
433
434	ENDDO
435
436	nqo = COUNT(delPhase(tracs(:)%name) == 'H2O' .AND. tracers_ico2lmdz(:)%component == 'lmdz') !--- Number of water phases
437
438	isotopes_ico2lmdz(1)%parent='H2O'
439	isotopes_ico2lmdz(1)%phase='gls'
440	isotopes_ico2lmdz(1)%nphas=3
441
442
443	nbtr=nqtot-nqo
444
445	ALLOCATE(conv_flg(nbtr))
446	ALLOCATE(pbl_flg(nbtr))
447
448	conv_flg(:) = 1 ! convection activated for all tracers
449	pbl_flg(:) = 1 ! boundary layer activated for all tracers
450
451
452	CALL init_infotrac_phy
453
454	! Initialize physical constant
455	day_length=86400
456	CALL getin('day_length',day_length)
457	CALL inifis(day_length,radius,g,kappa*cpp,cpp)
458
459
460
461	! init time
462	annee_ref=2015
463	CALL getin("anneeref",annee_ref)
464
465	day_ref=1
466	CALL getin("dayref",day_ref)
467
468	physics_timestep=dt*itau_physics
469	run_length=itaumax*dt
470	ndays=NINT(run_length/day_length)
471
472	day_ini=INT(itau0*dt/day_length)+day_ref
473	start_time= itau0dt/day_length-INT(itau0dt/day_length)
474
475	CALL init_time_lmdz(annee_ref, day_ref, day_ini, start_time, int(ndays), physics_timestep)
476
477	! Additional initializations for aquaplanets
478	! CALL getin("iflag_phys",iflag_phys)
479	! IF (iflag_phys>=100) THEN
480	! CALL iniaqua(klon_omp, iflag_phys)
481	! END IF
482
483
484
485	#ifdef INCA
486	CONTAINS
487
488	SUBROUTINE init_chem_trac()
489	IMPLICIT NONE
490
491	CALL Init_chem_inca_trac(nbtr)
492
493	END SUBROUTINE init_chem_trac
494
495	SUBROUTINE init_chem_transport()
496
497	IMPLICIT NONE
498
499	CALL init_transport(solsym, conv_flg,pbl_flg, hadv_inca, vadv_inca)
500
501	END SUBROUTINE init_chem_transport
502
503
504	#else
505	CONTAINS
506	SUBROUTINE init_chem_trac()
507	IMPLICIT NONE
508
509	END SUBROUTINE init_chem_trac
510
511	SUBROUTINE init_chem_transport()
512
513	IMPLICIT NONE
514
515	END SUBROUTINE init_chem_transport
516
517
518	#endif
519
520
521
522	END SUBROUTINE initialize_physics_omp
523
524
525
526
527	SUBROUTINE physics
528	USE icosa
529	USE time_mod
530	USE disvert_mod
531	USE transfert_mod
532	USE mpipara
533	USE omp_para
534	USE xios_mod
535	USE wxios
536	USE trace
537	USE distrib_icosa_lmdz_mod, ONLY : transfer_icosa_to_lmdz, transfer_lmdz_to_icosa
538	USE physics_external_mod, ONLY : it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q
539	USE write_field_mod
540	USE checksum_mod
541	USE vorticity_mod
542	USE tracer_icosa_mod
543	USE compute_transport_mod
544	USE, INTRINSIC :: IEEE_ARITHMETIC
545
546	! from LMDZ
547	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
548	USE geometry_mod, ONLY : cell_area
549	USE physiq_mod, ONLY: physiq
550	USE icolmdz_param_gravity_wave, ONLY: param_gravity_wave
551	USE mod_phys_lmdz_para, ONLY: reduce_sum, reduce_min
552	USE isotopes_mod, ONLY : iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO, toce, Rdefault, using_iso, initialisation_iso ! ==> LMDZ
553	USE write_field_phy
554	IMPLICIT NONE
555
556	REAL(rstd),POINTER :: phis(:)
557	REAL(rstd),POINTER :: ps(:)
558	REAL(rstd),POINTER :: theta_rhodz(:,:,:)
559	REAL(rstd),POINTER :: u(:,:)
560	REAL(rstd),POINTER :: wflux(:,:)
561	REAL(rstd),POINTER :: q(:,:,:)
562	REAL(rstd),POINTER :: p(:,:)
563	REAL(rstd),POINTER :: pks(:)
564	REAL(rstd),POINTER :: pk(:,:)
565	REAL(rstd),POINTER :: p_layer(:,:)
566	REAL(rstd),POINTER :: theta(:,:)
567	REAL(rstd),POINTER :: phi(:,:)
568	REAL(rstd),POINTER :: Temp(:,:)
569	REAL(rstd),POINTER :: ulon(:,:)
570	REAL(rstd),POINTER :: ulat(:,:)
571	REAL(rstd),POINTER :: vort(:,:)
572	REAL(rstd),POINTER :: vortc(:,:)
573	REAL(rstd),POINTER :: dulon(:,:)
574	REAL(rstd),POINTER :: dulat(:,:)
575	REAL(rstd),POINTER :: dTemp(:,:)
576	REAL(rstd),POINTER :: dq(:,:,:)
577	REAL(rstd),POINTER :: dps(:)
578	REAL(rstd),POINTER :: duc(:,:,:)
579
580
581	INTEGER :: ind,l,k, iq, iq_child
582
583	REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
584	!$OMP THREADPRIVATE(ps_phy)
585	REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
586	!$OMP THREADPRIVATE(p_phy)
587	REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
588	!$OMP THREADPRIVATE(p_layer_phy)
589	REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
590	!$OMP THREADPRIVATE(Temp_phy)
591	REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
592	!$OMP THREADPRIVATE(phis_phy)
593	REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
594	!$OMP THREADPRIVATE(phi_phy)
595	REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
596	!$OMP THREADPRIVATE(ulon_phy)
597	REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
598	!$OMP THREADPRIVATE(ulat_phy)
599	REAL(rstd),ALLOCATABLE,SAVE :: rot_phy(:,:)
600	!$OMP THREADPRIVATE(rot_phy)
601	REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
602	!$OMP THREADPRIVATE(q_phy)
603	REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
604	!$OMP THREADPRIVATE(wflux_phy)
605	REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
606	!$OMP THREADPRIVATE(dulon_phy)
607	REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
608	!$OMP THREADPRIVATE(dulat_phy)
609	REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
610	!$OMP THREADPRIVATE(dTemp_phy)
611	REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
612	!$OMP THREADPRIVATE(dq_phy)
613	REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
614	!$OMP THREADPRIVATE(dps_phy)
615	REAL(rstd) :: dtphy
616	LOGICAL :: debut
617	LOGICAL :: lafin
618	LOGICAL,SAVE :: first=.TRUE.
619	!$OMP THREADPRIVATE(first)
620	INTEGER, SAVE :: iter
621	!$OMP THREADPRIVATE(iter)
622	REAL(rstd),ALLOCATABLE,SAVE :: q_phy_old(:,:,:)
623	REAL :: qmin_omp,qmin_glo
624	REAL :: q_ave_bis(llm,nqtot)
625
626	IF(first) THEN
627	debut=.TRUE.
628	iter = 1
629	ELSE
630	debut=.FALSE.
631	iter=iter+1
632	ENDIF
633
634
635	IF(it-itau0>=itaumax) THEN
636	lafin=.TRUE.
637	ELSE
638	lafin=.FALSE.
639	ENDIF
640
641	IF (first) THEN
642	first=.FALSE.
643	CALL init_message(f_u,req_e1_vect,req_u)
644	CALL init_message(f_vort,req_z1_scal,req_z)
645	ALLOCATE(ps_phy(klon_omp))
646	ALLOCATE(p_phy(klon_omp,llm+1))
647	ALLOCATE(p_layer_phy(klon_omp,llm))
648	ALLOCATE(Temp_phy(klon_omp,llm))
649	ALLOCATE(phis_phy(klon_omp))
650	ALLOCATE(phi_phy(klon_omp,llm))
651	ALLOCATE(ulon_phy(klon_omp,llm))
652	ALLOCATE(ulat_phy(klon_omp,llm))
653	ALLOCATE(rot_phy(klon_omp,llm))
654	ALLOCATE(q_phy(klon_omp,llm,nqtot))
655	ALLOCATE(wflux_phy(klon_omp,llm))
656	ALLOCATE(dulon_phy(klon_omp,llm))
657	ALLOCATE(dulat_phy(klon_omp,llm))
658	ALLOCATE(dTemp_phy(klon_omp,llm))
659	ALLOCATE(dq_phy(klon_omp,llm,nqtot))
660	ALLOCATE(dps_phy(klon_omp))
661	ALLOCATE(q_phy_old(klon_omp,llm,nqtot))
662	!$OMP BARRIER
663
664	! CALL param_gravity_wave
665	ENDIF
666
667
668	!$OMP MASTER
669	! CALL update_calendar(it)
670	!$OMP END MASTER
671	!$OMP BARRIER
672	dtphy=itau_physics*dt
673
674
675
676	CALL transfert_message(f_u,req_u)
677	DO ind=1,ndomain
678	IF (assigned_domain(ind)) THEN
679	CALL swap_dimensions(ind)
680	CALL swap_geometry(ind)
681	u=f_u(ind)
682	vort=f_vort(ind)
683	CALL compute_vorticity(u,vort)
684	ENDIF
685	ENDDO
686
687	CALL transfert_message(f_vort,req_z)
688
689
690	DO ind=1,ndomain
691	CALL swap_dimensions(ind)
692	IF (assigned_domain(ind)) THEN
693	CALL swap_geometry(ind)
694
695	phis=f_phis(ind)
696	ps=f_ps(ind)
697	theta_rhodz=f_theta_rhodz(ind)
698	u=f_u(ind)
699	q=f_q(ind)
700	wflux=f_wflux(ind)
701	p=f_p(ind)
702	pks=f_pks(ind)
703	pk=f_pk(ind)
704	p_layer=f_p_layer(ind)
705	theta=f_theta(ind)
706	phi=f_phi(ind)
707	Temp=f_Temp(ind)
708	ulon=f_ulon(ind)
709	ulat=f_ulat(ind)
710	vort=f_vort(ind)
711	vortc=f_vortc(ind)
712
713	CALL grid_icosa_to_physics
714
715	ENDIF
716	ENDDO
717
718	!$OMP BARRIER
719	!$OMP MASTER
720	CALL SYSTEM_CLOCK(start_clock)
721	!$OMP END MASTER
722	CALL trace_start("physic")
723	! CALL trace_off()
724
725
726	! CALL writeField("p_in",f_p)
727	! CALL writeField("p_layer_in",f_p_layer)
728	! CALL writeField("phi_in",f_phi)
729	! CALL writeField("phis_in",f_phis)
730	! CALL writeField("ulon_in",f_ulon)
731	! CALL writeField("ulat_in",f_ulat)
732	! CALL writeField("Temp_in",f_Temp)
733	! CALL writeField("q_in",f_q)
734	! CALL writeField("wflux_in",f_wflux)
735	! CALL writeField("vortc",f_vortc)
736
737	! CALL checksum(f_p)
738	! CALL checksum(f_p_layer)
739	! CALL checksum(f_phi)
740	! CALL checksum(f_phis)
741	! CALL checksum(f_ulon)
742	! CALL checksum(f_ulat)
743	! CALL checksum(f_Temp)
744	! CALL checksum(f_q)
745	! CALL checksum(f_wflux)
746
747	CALL transfer_icosa_to_lmdz(f_p , p_phy)
748	CALL transfer_icosa_to_lmdz(f_p_layer, p_layer_phy)
749	CALL transfer_icosa_to_lmdz(f_phi , phi_phy)
750	CALL transfer_icosa_to_lmdz(f_phis , phis_phy )
751	CALL transfer_icosa_to_lmdz(f_ulon , ulon_phy )
752	CALL transfer_icosa_to_lmdz(f_ulat , ulat_phy)
753	CALL transfer_icosa_to_lmdz(f_vortc , rot_phy)
754	CALL transfer_icosa_to_lmdz(f_Temp , Temp_phy)
755	CALL transfer_icosa_to_lmdz(f_q , q_phy)
756	CALL transfer_icosa_to_lmdz(f_wflux , wflux_phy)
757
758	DO l=1,llm
759	wflux_phy(:,l) = - wflux_phy(:,l)*cell_area(:)
760	phi_phy(:,l)=phi_phy(:,l)-phis_phy(:)
761	ENDDO
762
763	DO iq=1, nqtot
764	IF (.NOT. tracers(iq)%has_parent) THEN
765	DO iq_child=1,tracers(iq)%nb_children
766	CALL tracers_dyn2phys(q_phy, tracers(iq)%children(iq_child), iq)
767	ENDDO
768	ENDIF
769
770	ENDDO
771
772	CALL wxios_set_context()
773
774	! Ehouarn: rot_phy() not implemented!! Set it to zero for now
775	! rot_phy(:,:)=0
776	CALL physiq(klon_omp, llm, debut, lafin, dtphy, &
777	p_phy, p_layer_phy, phi_phy, phis_phy, presnivs, &
778	ulon_phy, ulat_phy, rot_phy, Temp_phy, q_phy, wflux_phy, &
779	dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy)
780
781	CALL compute_q_ave(q_phy, q_ave)
782
783	DO iq=1, nqtot
784	q_phy(:,:,iq) = q_phy(:,:,iq) + dtphy*dq_phy(:,:,iq)
785	ENDDO
786
787	! IF (is_master) PRINT*, "=================================================="
788	! IF (is_master) PRINT*, "=================================================="
789	! CALL compute_q_ave(q_phy, q_ave_bis)
790	! DO iq=1, nqtot
791	! qmin_omp = MINVAL(q_phy(:,:,iq))
792	! CALL reduce_min(qmin_omp, qmin_glo)
793	! IF (is_master) PRINT*, "after 0 : MINVAL of ",TRIM(tracers(iq)%name), qmin_glo
794	! IF (is_master) PRINT*, "AVERAGE of ",TRIM(tracers(iq)%name), q_ave_bis(:,iq)
795	! ENDDO
796
797	CALL compute_qminimum(q_phy, q_ave, p_phy)
798
799	! DO iq=1, nqtot
800	! IF (.NOT. tracers(iq)%has_parent) THEN
801	! DO iq_child=1,tracers(iq)%nb_children
802	! CALL tracers_normalize(q_phy, tracers(iq)%children(iq_child), iq)
803	! ENDDO
804	! DO l=1,llm
805	! DO k=1,klon_omp
806	! IF ( q_phy(k,l,iq) < q_epsilonq_ave(l,iq)) q_phy(k,l,iq) = q_epsilonq_ave(l,iq)
807	! ENDDO
808	! ENDDO
809	! ENDIF
810	! ENDDO
811
812
813	! IF (is_master) PRINT*, "=================================================="
814	! IF (is_master) PRINT*, "=================================================="
815	!
816	! CALL compute_q_ave(q_phy, q_ave)
817	!
818	! DO iq=1, nqtot
819	! qmin_omp = MINVAL(q_phy(:,:,iq))
820	! CALL reduce_min(qmin_omp, qmin_glo)
821	! IF (is_master) PRINT*, "after 1 MINVAL of ",TRIM(tracers(iq)%name), qmin_glo
822	! IF (is_master) PRINT*, "AVERAGE of ",TRIM(tracers(iq)%name), q_ave(:,iq)
823	! IF (ANY(IEEE_IS_NAN(q_ave(:,iq)))) THEN
824	! CALL WriteField_phy("debug_"//TRIM(tracers(iq)%name),q_phy(:,:,iq),llm)
825	! IF (is_master) PRINT *,'NAN DETECTED'
826	! STOP
827	! ENDIF
828	! ENDDO
829
830	DO iq=1, nqtot
831	IF (.NOT. tracers(iq)%has_parent) THEN
832	DO iq_child=1,tracers(iq)%nb_children
833	CALL tracers_phys2dyn(q_phy, tracers(iq)%children(iq_child), iq)
834	ENDDO
835	ENDIF
836
837	ENDDO
838
839
840	CALL transfer_lmdz_to_icosa(dulon_phy, f_dulon )
841	CALL transfer_lmdz_to_icosa(dulat_phy, f_dulat )
842	CALL transfer_lmdz_to_icosa(dTemp_phy, f_dTemp )
843	CALL transfer_lmdz_to_icosa(q_phy , f_q )
844	CALL transfer_lmdz_to_icosa(dps_phy , f_dps )
845
846	! CALL writeField("dulon_out",f_dulon)
847	! CALL writeField("dulat_out",f_dulat)
848	! CALL writeField("dTemp_out",f_dTemp)
849	! CALL writeField("dq_out",f_dq)
850	! CALL writeField("dps_out",f_dps)
851
852	! CALL checksum(f_dulon)
853	! CALL checksum(f_dulat)
854	! CALL checksum(f_dTemp)
855	! CALL checksum(f_dq)
856	! CALL checksum(f_dps)
857
858	CALL send_message(f_dps,req_dps0)
859	CALL send_message(f_dulon,req_dulon0)
860	CALL send_message(f_dulat,req_dulat0)
861	CALL send_message(f_dTemp,req_dTemp0)
862	CALL send_message(f_q,req_dq0)
863
864	CALL wait_message(req_dps0)
865	CALL wait_message(req_dulon0)
866	CALL wait_message(req_dulat0)
867	CALL wait_message(req_dTemp0)
868	CALL wait_message(req_dq0)
869
870
871	! CALL trace_on()
872	CALL trace_end("physic")
873	!$OMP MASTER
874	CALL SYSTEM_CLOCK(stop_clock)
875	count_clock=count_clock+stop_clock-start_clock
876	!$OMP END MASTER
877
878	!$OMP BARRIER
879
880	DO ind=1,ndomain
881	CALL swap_dimensions(ind)
882	IF (assigned_domain(ind)) THEN
883	CALL swap_geometry(ind)
884
885	theta_rhodz=f_theta_rhodz(ind)
886	u=f_u(ind)
887	q=f_q(ind)
888	ps=f_ps(ind)
889	dulon=f_dulon(ind)
890	dulat=f_dulat(ind)
891	Temp=f_temp(ind)
892	dTemp=f_dTemp(ind)
893	dq=f_dq(ind)
894	dps=f_dps(ind)
895	duc=f_duc(ind)
896	p=f_p(ind)
897	pks=f_pks(ind)
898	pk=f_pk(ind)
899
900	CALL grid_physics_to_icosa
901	ENDIF
902	ENDDO
903
904	!$OMP BARRIER
905	CALL xios_set_context
906
907
908	CONTAINS
909
910	SUBROUTINE grid_icosa_to_physics
911	USE pression_mod
912	USE exner_mod
913	USE theta2theta_rhodz_mod
914	USE geopotential_mod
915	USE wind_from_lonlat_mod
916	USE omp_para
917	IMPLICIT NONE
918
919	REAL(rstd) :: uc(3)
920	INTEGER :: i,j,ij,l
921
922	! compute pression
923
924	DO l = ll_begin,ll_endp1
925	DO j=jj_begin,jj_end
926	DO i=ii_begin,ii_end
927	ij=(j-1)*iim+i
928	p(ij,l) = ap(l) + bp(l) * ps(ij)
929	ENDDO
930	ENDDO
931	ENDDO
932
933	!$OMP BARRIER
934
935	! compute exner
936
937	IF (is_omp_first_level) THEN
938	DO j=jj_begin,jj_end
939	DO i=ii_begin,ii_end
940	ij=(j-1)*iim+i
941	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
942	ENDDO
943	ENDDO
944	ENDIF
945
946	! 3D : pk
947	DO l = ll_begin,ll_end
948	DO j=jj_begin,jj_end
949	DO i=ii_begin,ii_end
950	ij=(j-1)*iim+i
951	pk(ij,l) = cpp * ((.5/preff)(p(ij,l)+p(ij,l+1))) * kappa
952	ENDDO
953	ENDDO
954	ENDDO
955
956	!$OMP BARRIER
957
958	! compute theta, temperature and pression at layer
959	DO l = ll_begin, ll_end
960	DO j=jj_begin,jj_end
961	DO i=ii_begin,ii_end
962	ij=(j-1)*iim+i
963	theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g)
964	Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp
965	p_layer(ij,l)=preff(pk(ij,l)/cpp)*(1./kappa)
966	ENDDO
967	ENDDO
968	ENDDO
969
970
971	!!! Compute geopotential
972
973	! for first layer
974	IF (is_omp_first_level) THEN
975	DO j=jj_begin,jj_end
976	DO i=ii_begin,ii_end
977	ij=(j-1)*iim+i
978	phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
979	ENDDO
980	ENDDO
981	ENDIF
982	!!-> implicit flush on phi(:,1)
983
984	!$OMP BARRIER
985
986	! for other layers
987	DO l = ll_beginp1, ll_end
988	DO j=jj_begin,jj_end
989	DO i=ii_begin,ii_end
990	ij=(j-1)*iim+i
991	phi(ij,l) = 0.5 * ( theta(ij,l) + theta(ij,l-1) ) &
992	* ( pk(ij,l-1) - pk(ij,l) )
993	ENDDO
994	ENDDO
995	ENDDO
996
997	!$OMP BARRIER
998
999
1000	IF (is_omp_first_level) THEN
1001	DO l = 2, llm
1002	DO j=jj_begin,jj_end
1003	! ---> Bug compilo intel ici en openmp
1004	! ---> Couper la boucle
1005	IF (j==jj_end+1) PRINT*,"this message must not be printed"
1006	DO i=ii_begin,ii_end
1007	ij=(j-1)*iim+i
1008	phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
1009	ENDDO
1010	ENDDO
1011	ENDDO
1012	! --> IMPLICIT FLUSH on phi --> non
1013	ENDIF
1014
1015	! compute wind centered lon lat compound
1016	DO l=ll_begin,ll_end
1017	DO j=jj_begin,jj_end
1018	DO i=ii_begin,ii_end
1019	ij=(j-1)*iim+i
1020	uc(:)=1/Ai(ij)* &
1021	( 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,:)) &
1022	+ 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,:)) &
1023	+ 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,:)) &
1024	+ 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,:)) &
1025	+ 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,:))&
1026	+ 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,:)))
1027	ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
1028	ulat(ij,l)=sum(uc(:)*elat_i(ij,:))
1029	ENDDO
1030	ENDDO
1031	ENDDO
1032
1033
1034	! compute centered vorticity
1035
1036	DO l=ll_begin,ll_end
1037	DO j=jj_begin,jj_end
1038	DO i=ii_begin,ii_end
1039	ij=(j-1)*iim+i
1040	vortc(ij,l) = Riv(ij,vup) * vort(ij+z_up,l) + &
1041	Riv(ij,vlup) * vort(ij+z_lup,l) + &
1042	Riv(ij,vldown)* vort(ij+z_ldown,l) + &
1043	Riv(ij,vdown) * vort(ij+z_down,l) + &
1044	Riv(ij,vrdown)* vort(ij+z_rdown,l) + &
1045	Riv(ij,vrup) * vort(ij+z_rup,l)
1046	ENDDO
1047	ENDDO
1048	ENDDO
1049
1050
1051
1052	!$OMP BARRIER
1053	END SUBROUTINE grid_icosa_to_physics
1054
1055
1056	SUBROUTINE grid_physics_to_icosa
1057	USE theta2theta_rhodz_mod
1058	USE omp_para
1059	IMPLICIT NONE
1060	INTEGER :: i,j,ij,l,iq
1061
1062	DO l=ll_begin,ll_end
1063	DO j=jj_begin,jj_end
1064	DO i=ii_begin,ii_end
1065	ij=(j-1)*iim+i
1066	duc(ij,:,l)=dulon(ij,l)elon_i(ij,:)+dulat(ij,l)elat_i(ij,:)
1067	ENDDO
1068	ENDDO
1069	ENDDO
1070
1071	DO l=ll_begin,ll_end
1072	DO j=jj_begin,jj_end
1073	DO i=ii_begin,ii_end
1074	ij=(j-1)*iim+i
1075	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,:) )
1076	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,:) )
1077	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,:) )
1078	ENDDO
1079	ENDDO
1080	ENDDO
1081
1082	DO l=ll_begin,ll_end
1083	DO j=jj_begin,jj_end
1084	DO i=ii_begin,ii_end
1085	ij=(j-1)*iim+i
1086	Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
1087	ENDDO
1088	ENDDO
1089	ENDDO
1090
1091	! DO iq=1,nqtot
1092	! DO l=ll_begin,ll_end
1093	! DO j=jj_begin,jj_end
1094	! DO i=ii_begin,ii_end
1095	! ij=(j-1)*iim+i
1096	! q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
1097	! ENDDO
1098	! ENDDO
1099	! ENDDO
1100	! ENDDO
1101
1102	!$OMP BARRIER
1103
1104	IF (is_omp_first_level) THEN
1105	DO j=jj_begin,jj_end
1106	DO i=ii_begin,ii_end
1107	ij=(j-1)*iim+i
1108	ps(ij)=ps(ij)+ dtphy * dps(ij)
1109	ENDDO
1110	ENDDO
1111	ENDIF
1112
1113	! CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
1114
1115	! compute pression
1116	!$OMP BARRIER
1117	DO l = ll_begin,ll_endp1
1118	DO j=jj_begin,jj_end
1119	DO i=ii_begin,ii_end
1120	ij=(j-1)*iim+i
1121	p(ij,l) = ap(l) + bp(l) * ps(ij)
1122	ENDDO
1123	ENDDO
1124	ENDDO
1125
1126	!$OMP BARRIER
1127
1128	! compute exner
1129
1130	IF (is_omp_first_level) THEN
1131	DO j=jj_begin,jj_end
1132	DO i=ii_begin,ii_end
1133	ij=(j-1)*iim+i
1134	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
1135	ENDDO
1136	ENDDO
1137	ENDIF
1138
1139	! 3D : pk
1140	DO l = ll_begin,ll_end
1141	DO j=jj_begin,jj_end
1142	DO i=ii_begin,ii_end
1143	ij=(j-1)*iim+i
1144	pk(ij,l) = cpp * ((.5/preff)(p(ij,l)+p(ij,l+1))) * kappa
1145	ENDDO
1146	ENDDO
1147	ENDDO
1148
1149	!$OMP BARRIER
1150
1151	! compute theta, temperature and pression at layer
1152	DO l = ll_begin, ll_end
1153	DO j=jj_begin,jj_end
1154	DO i=ii_begin,ii_end
1155	ij=(j-1)*iim+i
1156	theta_rhodz(ij,l,1) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
1157	ENDDO
1158	ENDDO
1159	ENDDO
1160
1161	END SUBROUTINE grid_physics_to_icosa
1162
1163	RECURSIVE SUBROUTINE tracers_dyn2phys(q, iq, iq_parent)
1164	USE tracer_icosa_mod
1165	IMPLICIT NONE
1166	REAL(rstd),INTENT(INOUT) :: q(:,:,:)
1167	INTEGER,INTENT(IN) :: iq, iq_parent
1168	INTEGER :: iq_child
1169
1170	q(:,:,iq) = q(:,:,iq)*q(:,:,iq_parent)
1171
1172	DO iq_child=1,tracers(iq)%nb_children
1173	CALL tracers_dyn2phys(q, tracers(iq)%children(iq_child), iq)
1174	ENDDO
1175
1176	END SUBROUTINE tracers_dyn2phys
1177
1178	RECURSIVE SUBROUTINE tracers_phys2dyn(q, iq, iq_parent)
1179	USE tracer_icosa_mod
1180	USE compute_transport_mod
1181	USE mod_phys_lmdz_para
1182	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
1183	USE dimphy
1184	USE mod_grid_phy_lmdz
1185	IMPLICIT NONE
1186	REAL(rstd),INTENT(INOUT) :: q(klon_omp,nbp_lev,nqtot)
1187	INTEGER,INTENT(IN) :: iq, iq_parent
1188	INTEGER :: iq_child
1189	INTEGER :: lij
1190
1191	DO iq_child=1,tracers(iq)%nb_children
1192	CALL tracers_phys2dyn(q, tracers(iq)%children(iq_child), iq)
1193	ENDDO
1194
1195	q(:,:,iq) = q(:,:,iq) / q(:,:,iq_parent)
1196
1197	END SUBROUTINE tracers_phys2dyn
1198
1199	RECURSIVE SUBROUTINE tracers_normalize(q, iq, iq_parent)
1200	USE tracer_icosa_mod
1201	USE compute_transport_mod
1202	USE mod_phys_lmdz_para
1203	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
1204	USE dimphy
1205	USE mod_grid_phy_lmdz
1206	IMPLICIT NONE
1207	REAL(rstd),INTENT(INOUT) :: q(klon_omp,nbp_lev,nqtot)
1208	INTEGER,INTENT(IN) :: iq, iq_parent
1209	INTEGER :: iq_child
1210	INTEGER :: l, k
1211
1212	DO iq_child=1,tracers(iq)%nb_children
1213	CALL tracers_normalize(q, tracers(iq)%children(iq_child), iq)
1214	ENDDO
1215
1216	DO l=1,klev
1217	DO k=1,klon
1218	IF ( q(k,l,iq_parent) >= q_epsilon*q_ave(l, iq_parent) ) THEN
1219	IF ( q(k,l,iq) < q_epsilon*q_ave(l, iq_parent) ) THEN
1220	q(k,l,iq) = q_ave(l, iq)/q_ave(l, iq_parent) * q(k,l,iq_parent)
1221	ENDIF
1222	ELSE
1223	q(k,l,iq) = q_ave(l, iq) * q_epsilon
1224	ENDIF
1225	ENDDO
1226	ENDDO
1227
1228	END SUBROUTINE tracers_normalize
1229
1230
1231	SUBROUTINE compute_q_ave(q, q_ave)
1232	USE dimphy
1233	USE tracer_icosa_mod
1234	USE mod_phys_lmdz_para
1235	USE mod_grid_phy_lmdz
1236	IMPLICIT NONE
1237	REAL(rstd),INTENT(IN) :: q(klon, klev ,nqtot)
1238	REAL(rstd),INTENT(OUT) :: q_ave(klev, nqtot)
1239	REAL(rstd) :: q_ave_omp(klev, nqtot)
1240	INTEGER :: iq, l , k
1241
1242	DO iq=1,nqtot
1243	DO l=1,klev
1244	q_ave_omp(l,iq) = 0.
1245	DO k=1,klon
1246	q_ave_omp(l,iq) = q_ave_omp(l,iq) + q(k,l,iq)
1247	ENDDO
1248	ENDDO
1249	ENDDO
1250	CALL reduce_sum(q_ave_omp, q_ave)
1251	CALL bcast(q_ave)
1252	q_ave(:,:) = q_ave(:,:) / klon_glo
1253	WHERE (q_ave < q_minimum) q_ave = q_minimum
1254
1255	END SUBROUTINE compute_q_ave
1256
1257	SUBROUTINE compute_qminimum(q, q_ave, p )
1258	USE dimphy
1259	USE tracer_icosa_mod
1260	IMPLICIT NONE
1261	REAL(rstd),INTENT(INOUT) :: q(klon, klev ,nqtot)
1262	REAL(rstd),INTENT(IN) :: q_ave(klev, nqtot)
1263	REAL(rstd),INTENT(IN) :: p(klon, klev)
1264	REAL(rstd) :: deltap(klon, klev)
1265	INTEGER :: iq, iq_g
1266	INTEGER :: iq_H2O_g, iq_H2O_l, iq_H216O_g, iq_H216O_l , iq_H218O_g, iq_H218O_l, iq_HDO_g, iq_HDO_l
1267
1268	DO iq=1, nqtot
1269	SELECT CASE (TRIM(tracers(iq)%name))
1270	CASE ("H2O_g") ; iq_H2O_g = iq
1271	CASE ("H2O_l") ; iq_H2O_l = iq
1272	CASE ("H216O_g") ; iq_H216O_g = iq
1273	CASE ("H216O_l") ; iq_H216O_l = iq
1274	CASE ("H218O_g") ; iq_H218O_g = iq
1275	CASE ("H218O_l") ; iq_H218O_l = iq
1276	CASE ("HDO_g") ; iq_HDO_g = iq
1277	CASE ("HDO_l") ; iq_HDO_l = iq
1278	CASE DEFAULT ; CYCLE
1279	END SELECT
1280	ENDDO
1281
1282	DO iq=1, nqtot
1283	IF (iq==iq_H2O_l) THEN
1284	iq_g = iq_H2O_g
1285	ELSE IF (iq==iq_H216O_l) THEN
1286	iq_g = iq_H216O_g
1287	ELSE IF (iq==iq_H218O_l) THEN
1288	iq_g = iq_H218O_g
1289	ELSE IF (iq==iq_HDO_l) THEN
1290	iq_g = iq_HDO_g
1291	ELSE
1292	CYCLE
1293	ENDIF
1294
1295	DO l=1,klev
1296	DO k=1,klon
1297	IF ( q(k,l,iq) <= q_epsilon*q_ave(l,iq)) THEN
1298	q(k,l,iq_g) = q(k,l,iq_g) - (q_epsilon*q_ave(l,iq) - q(k,l,iq))
1299	q(k,l,iq) = q_epsilon*q_ave(l,iq)
1300	ENDIF
1301	ENDDO
1302	ENDDO
1303	ENDDO
1304
1305	DO l=1,klev
1306	DO k=1,klon
1307	deltap(k,l) = p(k,l)-p(k,l+1)
1308	ENDDO
1309	ENDDO
1310
1311	DO iq=1, nqtot
1312	DO l=klev,2
1313	DO k=1,klon
1314	IF ( q(k,l,iq) < q_epsilon*q_ave(l,iq) ) THEN
1315	q(k,l-1,iq) = q(k,l-1,iq) - (q_epsilonq_ave(l,iq) - q(k,l,iq))deltap(k,l)/deltap(k,l-1)
1316	q(k,l,iq) = q_epsilon*q_ave(l,iq)
1317	ENDIF
1318	ENDDO
1319	ENDDO
1320	ENDDO
1321
1322	DO iq=1, nqtot
1323	DO l=1,klev-1
1324	DO k=1,klon
1325	IF ( q(k,l,iq) < q_epsilon*q_ave(l,iq) ) THEN
1326	q(k,l+1,iq) = q(k,l+1,iq) - (q_epsilonq_ave(l,iq) - q(k,l,iq))deltap(k,l)/deltap(k,l+1)
1327	q(k,l,iq) = q_epsilon*q_ave(l,iq)
1328	ENDIF
1329	ENDDO
1330
1331	DO k=1,klon
1332	IF (q(k,klev,iq) < q_epsilonq_ave(klev,iq)) q(k,klev,iq) = q_epsilonq_ave(klev,iq)
1333	ENDDO
1334	ENDDO
1335	ENDDO
1336
1337	END SUBROUTINE compute_qminimum
1338
1339
1340	END SUBROUTINE physics
1341
1342	END MODULE interface_icosa_lmdz_mod

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

Download in other formats:

Original Format