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: 27.7 KB

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

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

Download in other formats:

Original Format