Context Navigation

interface_icosa_lmdz.F90 @ 4536

Last change on this file since 4536 was 4501, checked in by Laurent Fairhead, 19 months ago

Modifications to the DYNAMICO / LMDZ interface needed for

the tracer parser
the new NetCDF96 library

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

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

Download in other formats:

Original Format