Context Navigation

interface_icosa_lmdz.f90 @ 2386

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

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

Download in other formats:

Original Format