Context Navigation

interface_icosa_lmdz.F90 @ 4498

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

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

Download in other formats:

Original Format