Context Navigation

interface_icosa_lmdz.f90 @ 3796

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

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

Download in other formats:

Original Format