Context Navigation

interface_icosa_lmdz.f90 @ 2005

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

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

Download in other formats:

Original Format