Context Navigation

interface_icosa_lmdz.f90 @ 3567

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

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

Download in other formats:

Original Format