Context Navigation

interface_icosa_lmdz.f90 @ 1862

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

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

Download in other formats:

Original Format