Context Navigation

interface_icosa_lmdz.f90 @ 3094

Last change on this file since 3094 was 2591, checked in by emillour, 3 years ago
Venus GCM: Add proper handling of "cpofT" case in the dynamico interface. EM
File size: 29.0 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	LOGICAL :: cpofT
202	REAL(rstd) :: nu_venus
203	REAL(rstd) :: t0_venus
204
205	! Calendar related stuff
206	REAL(rstd) :: ptimestep
207	REAL(rstd) :: run_length
208	INTEGER :: annee_ref
209	INTEGER :: day_ref
210	INTEGER :: day_ini
211	INTEGER :: day_end
212	INTEGER :: raz_date
213
214	! Tracer related information
215	INTEGER :: iflag_trac
216	CHARACTER(len=4) :: type_trac
217	! CHARACTER(len=30),ALLOCATABLE :: tname(:) ! tracer short name for restart and diagnostics
218	! CHARACTER(len=33),ALLOCATABLE :: ttext(:) ! tracer long name for diagnostics
219	TYPE(t_field),POINTER,SAVE :: f_ind_cell_glo(:)
220
221	INTEGER :: iflag_phys
222	INTEGER :: nq
223
224	!! to get values from startphy.nc controle array
225	!! ---------------------------------------------
226	logical :: startphy_file
227	! NetCDF stuff
228	integer :: status ! NetCDF return code
229	integer :: ncid ! NetCDF file ID
230	integer :: varid ! NetCDF variable ID
231	real :: tab_cntrl(100)
232
233	CALL init_distrib_icosa_lmdz
234
235	ALLOCATE(latfi(klon_omp))
236	ALLOCATE(lonfi(klon_omp))
237	ALLOCATE(airefi(klon_omp))
238	ALLOCATE(bounds_latfi(klon_omp,6))
239	ALLOCATE(bounds_lonfi(klon_omp,6))
240	ALLOCATE(ind_cell_glo(klon_omp))
241	ALLOCATE(dx(klon_omp))
242	ALLOCATE(dy(klon_omp))
243
244	CALL transfer_icosa_to_lmdz(geom%lat_i,latfi)
245	CALL transfer_icosa_to_lmdz(geom%lon_i,lonfi)
246	CALL transfer_icosa_to_lmdz(f_bounds_lat,bounds_latfi)
247	CALL transfer_icosa_to_lmdz(f_bounds_lon,bounds_lonfi)
248	CALL transfer_icosa_to_lmdz(geom%Ai,airefi)
249
250	CALL allocate_field(f_ind_cell_glo,field_t,type_real)
251
252	DO ind=1,ndomain
253	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
254	CALL swap_dimensions(ind)
255	CALL swap_geometry(ind)
256	DO j=jj_begin,jj_end
257	DO i=ii_begin,ii_end
258	ij=(j-1)*iim+i
259	f_ind_cell_glo(ind)%rval2d(ij)=domain(ind)%assign_cell_glo(i,j)
260	ENDDO
261	ENDDO
262	ENDDO
263
264
265	CALL transfer_icosa_to_lmdz(f_ind_cell_glo,ind_cell_glo)
266	CALL deallocate_field(f_ind_cell_glo)
267
268
269	! Initialize dimphy module
270	CALL init_dimphy(klon_omp,llm)
271
272	! Dummy initializations for dx(),dy() In principle these are not used
273	! in the physics; but this should be checked further...
274	dx(:)=1 ; dy(:)=1
275	CALL init_geometry(klon_omp,lonfi,latfi,bounds_lonfi,bounds_latfi,&
276	airefi,INT(ind_cell_glo),dx,dy)
277
278	scaleheight=scale_height/1000. ! Atmospheric scale height (km)
279	aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
280	bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
281	pseudoalt(:)=-scaleheight*log(presnivs(:)/preff)
282	CALL init_vertical_layers(llm,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
283
284	! Initialize planet_mod (quite redundant wrt vertical_levels...)
285	! CALL ini_planete_mod(llm,preff,ap,bp)
286
287	! Initialize tracer names, numbers, etc. for physics
288
289	! init tracers model for standard lmdz case
290	!$OMP MASTER
291	ALLOCATE(tname(nqtot))
292	ALLOCATE(ttext(nqtot))
293	!$OMP END MASTER
294	!$OMP BARRIER
295
296	! read tname() from traceur.def file
297	IF (is_mpi_root) THEN
298	!$OMP MASTER
299	OPEN(unit=42,file="traceur.def",form="formatted",status="old",iostat=ierr)
300	IF (ierr==0) THEN
301	READ(42,*) nq ! should be the same as nqtot
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	DO i=1,nqtot
307	READ(42,*) j,k,tname(i)
308	ttext(i)=trim(tname(i))//"VL1"
309	ENDDO
310	CLOSE(42)
311	ENDIF
312	!$OMP END MASTER
313	!$OMP BARRIER
314	ENDIF ! of (is_mpi_root)
315
316	DO i=1,nqtot
317	CALL bcast(tname(i))
318	CALL bcast(ttext(i))
319	ENDDO
320
321	! Get/set some constants for the physics
322
323
324	startphy_file=.true.
325	CALL getin('startphy_file',startphy_file)
326
327	! value in physics daysec=10087066.76s, important for solar radiation,
328	! for the rest ptime/timeofday is what matters. To well estimate them,
329	! day_length must be multiple of dt*itau_physics. Error order:1e-5.
330	day_length=10087000
331	CALL getin('day_length',day_length)
332
333	run_length=day_length ! default
334	CALL getin('run_length',run_length)
335
336	raz_date=0 ! default: 0: no change in date
337	CALL getin('raz_date',raz_date)
338
339	IF (startphy_file) THEN
340	! Read in some information from the startphy.nc file
341
342	IF (is_mpi_root) THEN
343	!$OMP MASTER
344	status=nf90_open('startphy.nc',NF90_NOWRITE,ncid)
345	if (status.ne.nf90_noerr) then
346	write(,)"Failed to open startphy.nc"
347	write(,)trim(nf90_strerror(status))
348	stop
349	endif
350
351	status=nf90_inq_varid(ncid,"controle",varid)
352	if (status.ne.nf90_noerr) then
353	write(,)"Failed to find controle variable"
354	write(,)trim(nf90_strerror(status))
355	stop
356	endif
357
358	status=nf90_get_var(ncid,varid,tab_cntrl)
359	! extract needed variables from tab_cntrl
360	day_ini=tab_cntrl(13)
361	annee_ref=tab_cntrl(14)
362	! check if tab_cntrl(1), the stored physics time step
363	! is the same as the current physics time step (within roundoff precision)
364	if (abs(((itau_physicsdt)-tab_cntrl(1))/(itau_physicsdt))<=1.e-8) then
365	! Everything OK
366	ptime=modulo((tab_cntrl(15)*tab_cntrl(1))/day_length,1.0)
367	else ! unless raz_date == 1 , we have a problem
368	if (raz_date==1) then
369	! we reset date to midnight at lon=0 in the physics
370	ptime=0.0
371	day_ini=0
372	else
373	write(,)"Error: physics time step in startphy.nc is different"
374	write(,)" from what is specified via run_icosa.def"
375	write(,)" From run_icosa.def: ",itau_physics*dt
376	write(,)" From startphy.nc file: ",tab_cntrl(1)
377	write(,)" You must reset date to midnight at lon=0"
378	write(,)" by specifying raz_date=1 in your def file"
379	call abort
380	endif
381	endif
382
383	status=nf90_close(ncid)
384	!$OMP END MASTER
385	!$OMP BARRIER
386	ENDIF ! of !IF (is_mpi_root)
387
388	CALL bcast(day_ini)
389	CALL bcast(annee_ref)
390	CALL bcast(ptime)
391
392	ELSE
393	! required information that is not in tab_cntrl
394	! has to be default or read from def files
395	day_ini=0
396	annee_ref=1
397	CALL getin('annee_ref',annee_ref)
398	ptime=0.
399	ENDIF
400
401	day_end=day_ini+nint(run_length/day_length)
402
403	! Other required values which have to be read from def files
404	day_ref=1
405	CALL getin('day_ref',day_ref)
406	iflag_trac=0
407	CALL getin('iflag_trac',iflag_trac)
408
409
410	! Initialize some physical constants
411	CALL suphec(cpp)
412
413	! Initialize cpdet_phy module
414	! FLAG CP CONSTANT OR NOT... (in LMDZ.COMMON: cpofT in gcm.def)
415	! here, it allows to run the physics with a variable Cp together with dynamics
416	! with constant Cp.
417	! Conversion Theta/T takes this into account (see below)
418	cpofT=.false. ! default, assume Cp to be constant in physics
419	CALL getin('cpofT',cpofT)
420
421	IF (cpofT) THEN
422	nu_venus=0.35
423	t0_venus=460.
424	ELSE
425	nu_venus=0.
426	t0_venus=0.
427	ENDIF
428
429	CALL init_cpdet_phy(cpp,nu_venus,t0_venus)
430
431	! Initialize some "temporal and calendar" related variables
432	ptimestep=itau_physics*dt
433	CALL init_time_lmdz(annee_ref,day_ref,day_ini,day_end,ptimestep)
434
435	! Initialize tracers in physics
436	CALL init_infotrac_phy(iflag_trac,nqtot,tname,ttext)
437
438
439	! Initializations of some module variables
440
441	!$OMP MASTER
442	! initialize pday
443	pday = day_ini
444	!$OMP END MASTER
445	!$OMP BARRIER
446
447	END SUBROUTINE initialize_physics_omp
448
449
450
451
452	SUBROUTINE physics
453	USE icosa
454	USE time_mod
455	USE disvert_mod
456	USE transfert_mod
457	USE mpipara
458	USE xios_mod
459	USE wxios
460	USE trace
461	USE distrib_icosa_lmdz_mod, ONLY : transfer_icosa_to_lmdz, transfer_lmdz_to_icosa
462	USE physics_external_mod, ONLY : it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q
463	USE write_field_mod
464	USE checksum_mod
465	! from LMDZ
466	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
467	USE geometry_mod, ONLY : cell_area
468	USE physiq_mod, ONLY: physiq
469	IMPLICIT NONE
470
471	REAL(rstd),POINTER :: phis(:)
472	REAL(rstd),POINTER :: ps(:)
473	REAL(rstd),POINTER :: theta_rhodz(:,:,:)
474	REAL(rstd),POINTER :: u(:,:)
475	REAL(rstd),POINTER :: wflux(:,:)
476	REAL(rstd),POINTER :: q(:,:,:)
477	REAL(rstd),POINTER :: p(:,:)
478	REAL(rstd),POINTER :: pks(:)
479	REAL(rstd),POINTER :: pk(:,:)
480	REAL(rstd),POINTER :: p_layer(:,:)
481	REAL(rstd),POINTER :: theta(:,:)
482	REAL(rstd),POINTER :: phi(:,:)
483	REAL(rstd),POINTER :: Temp(:,:)
484	REAL(rstd),POINTER :: ulon(:,:)
485	REAL(rstd),POINTER :: ulat(:,:)
486	REAL(rstd),POINTER :: dulon(:,:)
487	REAL(rstd),POINTER :: dulat(:,:)
488	REAL(rstd),POINTER :: dTemp(:,:)
489	REAL(rstd),POINTER :: dq(:,:,:)
490	REAL(rstd),POINTER :: dps(:)
491	REAL(rstd),POINTER :: duc(:,:,:)
492
493
494	INTEGER :: ind,l
495
496	REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
497	!$OMP THREADPRIVATE(ps_phy)
498	REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
499	!$OMP THREADPRIVATE(p_phy)
500	REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
501	!$OMP THREADPRIVATE(p_layer_phy)
502	REAL(rstd),ALLOCATABLE,SAVE :: pk_phy(:,:)
503	!$OMP THREADPRIVATE(pk_phy)
504	REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
505	!$OMP THREADPRIVATE(Temp_phy)
506	REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
507	!$OMP THREADPRIVATE(phis_phy)
508	REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
509	!$OMP THREADPRIVATE(phi_phy)
510	REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
511	!$OMP THREADPRIVATE(ulon_phy)
512	REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
513	!$OMP THREADPRIVATE(ulat_phy)
514	REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
515	!$OMP THREADPRIVATE(q_phy)
516	REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
517	!$OMP THREADPRIVATE(wflux_phy)
518	REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
519	!$OMP THREADPRIVATE(dulon_phy)
520	REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
521	!$OMP THREADPRIVATE(dulat_phy)
522	REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
523	!$OMP THREADPRIVATE(dTemp_phy)
524	REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
525	!$OMP THREADPRIVATE(dq_phy)
526	REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
527	!$OMP THREADPRIVATE(dps_phy)
528	REAL(rstd) :: dtphy
529	LOGICAL :: debut
530	LOGICAL :: lafin
531	LOGICAL,SAVE :: first=.TRUE.
532	!$OMP THREADPRIVATE(first)
533
534	IF(first) THEN
535	debut=.TRUE.
536	ELSE
537	debut=.FALSE.
538	ENDIF
539
540
541	IF(it-itau0>=itaumax) THEN
542	lafin=.TRUE.
543	ELSE
544	lafin=.FALSE.
545	ENDIF
546
547	IF (first) THEN
548	first=.FALSE.
549	CALL init_message(f_u,req_e1_vect,req_u)
550	ALLOCATE(ps_phy(klon_omp))
551	ALLOCATE(p_phy(klon_omp,llm+1))
552	ALLOCATE(p_layer_phy(klon_omp,llm))
553	ALLOCATE(pk_phy(klon_omp,llm))
554	ALLOCATE(Temp_phy(klon_omp,llm))
555	ALLOCATE(phis_phy(klon_omp))
556	ALLOCATE(phi_phy(klon_omp,llm))
557	ALLOCATE(ulon_phy(klon_omp,llm))
558	ALLOCATE(ulat_phy(klon_omp,llm))
559	ALLOCATE(q_phy(klon_omp,llm,nqtot))
560	ALLOCATE(wflux_phy(klon_omp,llm))
561	ALLOCATE(dulon_phy(klon_omp,llm))
562	ALLOCATE(dulat_phy(klon_omp,llm))
563	ALLOCATE(dTemp_phy(klon_omp,llm))
564	ALLOCATE(dq_phy(klon_omp,llm,nqtot))
565	ALLOCATE(dps_phy(klon_omp))
566	!$OMP BARRIER
567	ENDIF
568
569
570	!$OMP MASTER
571	! CALL update_calendar(it)
572	!$OMP END MASTER
573	!$OMP BARRIER
574	dtphy=itau_physics*dt
575
576
577
578	CALL transfert_message(f_u,req_u)
579
580	DO ind=1,ndomain
581	CALL swap_dimensions(ind)
582	IF (assigned_domain(ind)) THEN
583	CALL swap_geometry(ind)
584
585	phis=f_phis(ind)
586	ps=f_ps(ind)
587	theta_rhodz=f_theta_rhodz(ind)
588	u=f_u(ind)
589	q=f_q(ind)
590	wflux=f_wflux(ind)
591	p=f_p(ind)
592	pks=f_pks(ind)
593	pk=f_pk(ind)
594	p_layer=f_p_layer(ind)
595	theta=f_theta(ind)
596	phi=f_phi(ind)
597	Temp=f_Temp(ind)
598	ulon=f_ulon(ind)
599	ulat=f_ulat(ind)
600
601	CALL grid_icosa_to_physics
602
603	ENDIF
604	ENDDO
605
606	!$OMP BARRIER
607	!$OMP MASTER
608	CALL SYSTEM_CLOCK(start_clock)
609	!$OMP END MASTER
610	CALL trace_start("physic")
611	! CALL trace_off()
612
613
614	! CALL writeField("p_in",f_p)
615	! CALL writeField("p_layer_in",f_p_layer)
616	! CALL writeField("phi_in",f_phi)
617	! CALL writeField("phis_in",f_phis)
618	! CALL writeField("ulon_in",f_ulon)
619	! CALL writeField("ulat_in",f_ulat)
620	! CALL writeField("Temp_in",f_Temp)
621	! CALL writeField("q_in",f_q)
622	! CALL writeField("wflux_in",f_wflux)
623
624	! CALL checksum(f_p)
625	! CALL checksum(f_p_layer)
626	! CALL checksum(f_phi)
627	! CALL checksum(f_phis)
628	! CALL checksum(f_ulon)
629	! CALL checksum(f_ulat)
630	! CALL checksum(f_Temp)
631	! CALL checksum(f_q)
632	! CALL checksum(f_wflux)
633
634	CALL transfer_icosa_to_lmdz(f_p , p_phy)
635	CALL transfer_icosa_to_lmdz(f_p_layer, p_layer_phy)
636	CALL transfer_icosa_to_lmdz(f_pk , pk_phy)
637	CALL transfer_icosa_to_lmdz(f_phi , phi_phy)
638	CALL transfer_icosa_to_lmdz(f_phis , phis_phy )
639	CALL transfer_icosa_to_lmdz(f_ulon , ulon_phy )
640	CALL transfer_icosa_to_lmdz(f_ulat , ulat_phy)
641	CALL transfer_icosa_to_lmdz(f_Temp , Temp_phy)
642	CALL transfer_icosa_to_lmdz(f_q , q_phy)
643	CALL transfer_icosa_to_lmdz(f_wflux , wflux_phy)
644
645	DO l=1,llm
646	! Warning: In the physics, vertical flux convention is positive if downwards!
647	wflux_phy(:,l)= - wflux_phy(:,l)*cell_area(:)
648	! Compute relative geopotential
649	phi_phy(:,l)=phi_phy(:,l)-phis_phy(:)
650	ENDDO
651
652	CALL wxios_set_context()
653
654	! Update pday and ptime to send to physics
655	!$OMP MASTER
656	ptime=ptime+dtphy/day_length
657	IF (ptime >= 1.) THEN
658	ptime=ptime-1
659	pday=pday+1
660	ENDIF
661	!$OMP END MASTER
662	!$OMP BARRIER
663
664	CALL physiq(klon_omp, llm, nqtot, &
665	debut, lafin, &
666	pday, ptime, dtphy, &
667	p_phy, p_layer_phy, pk_phy, &
668	phi_phy, phis_phy, presnivs, &
669	ulon_phy, ulat_phy, Temp_phy, q_phy, &
670	wflux_phy, &
671	dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy)
672
673	CALL transfer_lmdz_to_icosa(dulon_phy, f_dulon )
674	CALL transfer_lmdz_to_icosa(dulat_phy, f_dulat )
675	CALL transfer_lmdz_to_icosa(dTemp_phy, f_dTemp )
676	CALL transfer_lmdz_to_icosa(dq_phy , f_dq )
677	CALL transfer_lmdz_to_icosa(dps_phy , f_dps )
678
679	! CALL writeField("dulon_out",f_dulon)
680	! CALL writeField("dulat_out",f_dulat)
681	! CALL writeField("dTemp_out",f_dTemp)
682	! CALL writeField("dq_out",f_dq)
683	! CALL writeField("dps_out",f_dps)
684
685	! CALL checksum(f_dulon)
686	! CALL checksum(f_dulat)
687	! CALL checksum(f_dTemp)
688	! CALL checksum(f_dq)
689	! CALL checksum(f_dps)
690
691	CALL send_message(f_dps,req_dps0)
692	CALL send_message(f_dulon,req_dulon0)
693	CALL send_message(f_dulat,req_dulat0)
694	CALL send_message(f_dTemp,req_dTemp0)
695	CALL send_message(f_dq,req_dq0)
696
697	CALL wait_message(req_dps0)
698	CALL wait_message(req_dulon0)
699	CALL wait_message(req_dulat0)
700	CALL wait_message(req_dTemp0)
701	CALL wait_message(req_dq0)
702
703
704	! CALL trace_on()
705	CALL trace_end("physic")
706	!$OMP MASTER
707	CALL SYSTEM_CLOCK(stop_clock)
708	count_clock=count_clock+stop_clock-start_clock
709	!$OMP END MASTER
710
711	!$OMP BARRIER
712
713	DO ind=1,ndomain
714	CALL swap_dimensions(ind)
715	IF (assigned_domain(ind)) THEN
716	CALL swap_geometry(ind)
717
718	theta_rhodz=f_theta_rhodz(ind)
719	u=f_u(ind)
720	q=f_q(ind)
721	ps=f_ps(ind)
722	dulon=f_dulon(ind)
723	dulat=f_dulat(ind)
724	Temp=f_temp(ind)
725	dTemp=f_dTemp(ind)
726	dq=f_dq(ind)
727	dps=f_dps(ind)
728	duc=f_duc(ind)
729	p=f_p(ind)
730	pks=f_pks(ind)
731	pk=f_pk(ind)
732
733	CALL grid_physics_to_icosa
734	ENDIF
735	ENDDO
736
737	!$OMP BARRIER
738	CALL xios_set_context
739
740
741	CONTAINS
742
743	SUBROUTINE grid_icosa_to_physics
744	USE pression_mod
745	USE exner_mod
746	USE theta2theta_rhodz_mod
747	USE geopotential_mod
748	USE wind_mod
749	USE omp_para
750	USE cpdet_phy_mod, only: tpot2t
751	IMPLICIT NONE
752
753	REAL(rstd) :: uc(3)
754	INTEGER :: i,j,ij,l
755	REAL :: tmpT(1),tmptheta(1),tmppk(1)
756
757
758	! compute pression
759
760	DO l = ll_begin,ll_endp1
761	DO j=jj_begin,jj_end
762	DO i=ii_begin,ii_end
763	ij=(j-1)*iim+i
764	p(ij,l) = ap(l) + bp(l) * ps(ij)
765	ENDDO
766	ENDDO
767	ENDDO
768
769	!$OMP BARRIER
770
771	! compute exner
772
773	IF (is_omp_first_level) THEN
774	DO j=jj_begin,jj_end
775	DO i=ii_begin,ii_end
776	ij=(j-1)*iim+i
777	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
778	ENDDO
779	ENDDO
780	ENDIF
781
782	! 3D : pk
783	DO l = ll_begin,ll_end
784	DO j=jj_begin,jj_end
785	DO i=ii_begin,ii_end
786	ij=(j-1)*iim+i
787	pk(ij,l) = cpp * ((.5/preff)(p(ij,l)+p(ij,l+1))) * kappa
788	ENDDO
789	ENDDO
790	ENDDO
791
792	!$OMP BARRIER
793
794	! compute theta, temperature and pression at layer
795	! Conversion Theta/T takes Cp(T) into account
796	DO l = ll_begin, ll_end
797	DO j=jj_begin,jj_end
798	DO i=ii_begin,ii_end
799	ij=(j-1)*iim+i
800	theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g)
801	tmptheta(1) = theta(ij,l)
802	tmppk(1) = pk(ij,l)
803	call tpot2t(1,tmptheta,tmpT,tmppk)
804	Temp(ij,l) = tmpT(1)
805	p_layer(ij,l)=preff(pk(ij,l)/cpp)*(1./kappa)
806	ENDDO
807	ENDDO
808	ENDDO
809
810
811	!!! Compute geopotential
812	!! This computation (with teta = cp T / pk !) is identical to
813	!! delta phi = R/RMD T/p delta p (r=R/RMD=cpp*kappa)
814
815	! for first layer
816	IF (is_omp_first_level) THEN
817	DO j=jj_begin,jj_end
818	DO i=ii_begin,ii_end
819	ij=(j-1)*iim+i
820	phi( ij,1 ) = phis( ij ) + ( cpp*Temp(ij,1)/pk(ij,1) ) &
821	* ( pks(ij) - pk(ij,1) )
822	ENDDO
823	ENDDO
824	ENDIF
825	!!-> implicit flush on phi(:,1)
826
827	! for other layers
828	DO l = ll_beginp1, ll_end
829	DO j=jj_begin,jj_end
830	DO i=ii_begin,ii_end
831	ij=(j-1)*iim+i
832	phi(ij,l) = 0.5 cpp( Temp(ij,l)/pk(ij,l)+Temp(ij,l-1)/pk(ij,l-1) ) &
833	* ( pk(ij,l-1) - pk(ij,l) )
834	ENDDO
835	ENDDO
836	ENDDO
837
838	!$OMP BARRIER
839
840
841	IF (is_omp_first_level) THEN
842	DO l = 2, llm
843	DO j=jj_begin,jj_end
844	! ---> Bug compilo intel ici en openmp
845	! ---> Couper la boucle
846	IF (j==jj_end+1) PRINT*,"this message must not be printed"
847	DO i=ii_begin,ii_end
848	ij=(j-1)*iim+i
849	phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
850	ENDDO
851	ENDDO
852	ENDDO
853	! --> IMPLICIT FLUSH on phi --> non
854	ENDIF
855
856	! compute wind centered lon lat compound
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	uc(:)=1/Ai(ij)* &
862	( 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,:)) &
863	+ 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,:)) &
864	+ 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,:)) &
865	+ 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,:)) &
866	+ 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,:))&
867	+ 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,:)))
868	ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
869	ulat(ij,l)=sum(uc(:)*elat_i(ij,:))
870	ENDDO
871	ENDDO
872	ENDDO
873
874	!$OMP BARRIER
875	END SUBROUTINE grid_icosa_to_physics
876
877
878	SUBROUTINE grid_physics_to_icosa
879	USE theta2theta_rhodz_mod
880	USE omp_para
881	USE cpdet_phy_mod, only: t2tpot
882	IMPLICIT NONE
883	INTEGER :: i,j,ij,l,iq
884	REAL :: tmptheta(1),tmpT(1),tmppk(1)
885
886	DO l=ll_begin,ll_end
887	DO j=jj_begin,jj_end
888	DO i=ii_begin,ii_end
889	ij=(j-1)*iim+i
890	duc(ij,:,l)=dulon(ij,l)elon_i(ij,:)+dulat(ij,l)elat_i(ij,:)
891	ENDDO
892	ENDDO
893	ENDDO
894
895	DO l=ll_begin,ll_end
896	DO j=jj_begin,jj_end
897	DO i=ii_begin,ii_end
898	ij=(j-1)*iim+i
899	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,:) )
900	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,:) )
901	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,:) )
902	ENDDO
903	ENDDO
904	ENDDO
905
906	DO l=ll_begin,ll_end
907	DO j=jj_begin,jj_end
908	DO i=ii_begin,ii_end
909	ij=(j-1)*iim+i
910	Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
911	ENDDO
912	ENDDO
913	ENDDO
914
915	DO iq=1,nqtot
916	DO l=ll_begin,ll_end
917	DO j=jj_begin,jj_end
918	DO i=ii_begin,ii_end
919	ij=(j-1)*iim+i
920	q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
921	ENDDO
922	ENDDO
923	ENDDO
924	ENDDO
925
926	!$OMP BARRIER
927
928	IF (is_omp_first_level) THEN
929	DO j=jj_begin,jj_end
930	DO i=ii_begin,ii_end
931	ij=(j-1)*iim+i
932	ps(ij)=ps(ij)+ dtphy * dps(ij)
933	ENDDO
934	ENDDO
935	ENDIF
936
937	! CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
938
939	! compute pression
940	!$OMP BARRIER
941	DO l = ll_begin,ll_endp1
942	DO j=jj_begin,jj_end
943	DO i=ii_begin,ii_end
944	ij=(j-1)*iim+i
945	p(ij,l) = ap(l) + bp(l) * ps(ij)
946	ENDDO
947	ENDDO
948	ENDDO
949
950	!$OMP BARRIER
951
952	! compute exner
953
954	IF (is_omp_first_level) THEN
955	DO j=jj_begin,jj_end
956	DO i=ii_begin,ii_end
957	ij=(j-1)*iim+i
958	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
959	ENDDO
960	ENDDO
961	ENDIF
962
963	! 3D : pk
964	DO l = ll_begin,ll_end
965	DO j=jj_begin,jj_end
966	DO i=ii_begin,ii_end
967	ij=(j-1)*iim+i
968	pk(ij,l) = cpp * ((.5/preff)(p(ij,l)+p(ij,l+1))) * kappa
969	ENDDO
970	ENDDO
971	ENDDO
972
973	!$OMP BARRIER
974
975	! compute theta, temperature and pression at layer
976	! Conversion Theta/T takes Cp(T) into account
977	DO l = ll_begin, ll_end
978	DO j=jj_begin,jj_end
979	DO i=ii_begin,ii_end
980	ij=(j-1)*iim+i
981	tmpT(1) = Temp(ij,l)
982	tmppk(1) = pk(ij,l)
983	call t2tpot(1,tmpT,tmptheta,tmppk)
984	theta_rhodz(ij,l,1) = tmptheta(1) * ((p(ij,l)-p(ij,l+1))/g)
985	ENDDO
986	ENDDO
987	ENDDO
988
989	END SUBROUTINE grid_physics_to_icosa
990
991
992
993	END SUBROUTINE physics
994
995
996
997
998
999	END MODULE interface_icosa_lmdz_mod

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

Download in other formats:

Original Format