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

Rev	Line
[1696]	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
[1710]	174
	175	USE netcdf
[1696]	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
[2785]	206
[1696]	207	TYPE(t_field),POINTER,SAVE :: f_ind_cell_glo(:)
	208
	209	INTEGER :: iflag_phys
	210	INTEGER :: nq
[2785]	211	character(len=500) :: tracline ! to read traceur.def file
[1696]	212
[1710]	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
[1696]	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))
[1698]	264	pseudoalt(:)=-scaleheight*log(presnivs(:)/preff)
[1696]	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
[2785]	282	! read nq from traceur.def file
	283	! We allow to read the Modern-trac-V1 as well to get nqtot
[2000]	284	IF (is_mpi_root) THEN
	285	!$OMP MASTER
[1696]	286	OPEN(unit=42,file="traceur.def",form="formatted",status="old",iostat=ierr)
	287	IF (ierr==0) THEN
[2785]	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')
[1696]	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)
[2785]	307	ENDIF !ierr==0 (1st one)
[1696]	308	!$OMP END MASTER
	309	!$OMP BARRIER
[2000]	310	ENDIF ! of (is_mpi_root)
[1696]	311
[2000]	312
[2785]	313
[1696]	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)
[1710]	325
	326	startphy_file=.true.
	327	CALL getin('startphy_file',startphy_file)
	328
	329	IF (startphy_file) THEN
	330
[2000]	331	IF (is_mpi_root) THEN
	332	!$OMP MASTER
[1710]	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)
[2000]	349	print*,"initialize_physics_omp: start_day",start_day
[1710]	350	status=nf90_close(ncid)
[2000]	351	!$OMP END MASTER
	352	!$OMP BARRIER
	353	ENDIF ! of !IF (is_mpi_root)
	354
	355	CALL bcast(start_day)
[1710]	356
	357	ELSE
	358
	359	start_day=0
	360	CALL getin('start_day',start_day)
	361
	362	ENDIF
	363
[1696]	364	day_length=86400
	365	CALL getin('day_length',day_length)
[1710]	366	year_length = 31557600 ! 365.25 * 86400
	367	CALL getin('year_length',year_length)
[1696]	368	ndays=nint(itaumax*(dt/day_length))! number of days to run
	369	physics_timestep=dt*itau_physics
[2000]	370
[2556]	371	CALL inifis(klon_omp,llm,nqtot,start_day,day_length,nint(ndays),physics_timestep, &
[1696]	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
[2000]	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
[1696]	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