Context Navigation

interface_icosa_lmdz.f90 @ 2368

Last change on this file since 2368 was 2000, checked in by emillour, 7 years ago

Dynamico-physics interfaces:

some OpenMP fixes.
small bug fix: in Dynamico vertical flux is positive when upwards, but in the physics positive is downwards, so change the sign when sending values to the physics.

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

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

Download in other formats:

Original Format