Context Navigation

interface_icosa_lmdz.f90 @ 3094

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

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

Download in other formats:

Original Format