Context Navigation

interface_icosa_lmdz.f90 @ 1847

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

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

Download in other formats:

Original Format