Context Navigation

interface_icosa_lmdz.F90 @ 5163

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

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

Download in other formats:

Original Format