Context Navigation

interface_icosa_lmdz.F90 @ 5163

Last change on this file since 5163 was 5059, checked in by idelkadi, 5 months ago
Correction to correct: ambiguity for the 2 subroutines etat0 and init_etat20, which are defined in 2 different modules.
File size: 32.1 KB

Rev	Line
[4855]	1	MODULE interface_icosa_lmdz_mod
[4498]	2
	3	USE field_mod, ONLY: t_field
	4	USE transfert_mod, ONLY: t_message
	5
	6
	7	TYPE(t_message),SAVE :: req_u, req_z
	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_vort(:)
	20	TYPE(t_field),POINTER,SAVE :: f_vortc(:)
	21	TYPE(t_field),POINTER,SAVE :: f_dulon(:)
	22	TYPE(t_field),POINTER,SAVE :: f_dulat(:)
	23	TYPE(t_field),POINTER,SAVE :: f_dTemp(:)
	24	TYPE(t_field),POINTER,SAVE :: f_dq(:)
	25	TYPE(t_field),POINTER,SAVE :: f_dps(:)
	26	TYPE(t_field),POINTER,SAVE :: f_duc(:)
	27	TYPE(t_field),POINTER,SAVE :: f_bounds_lon(:)
	28	TYPE(t_field),POINTER,SAVE :: f_bounds_lat(:)
	29
	30	INTEGER :: start_clock
	31	INTEGER :: stop_clock
	32	INTEGER :: count_clock=0
	33
	34	INTEGER,SAVE :: nbp_phys
	35	INTEGER,SAVE :: nbp_phys_glo
	36
	37
	38	CONTAINS
	39
[4855]	40	SUBROUTINE pre_initialize_physics
[5059]	41	USE etat0_plugin_mod, ONLY : init_etat0_plugin, etat0_plugin
	42	USE icolmdz_etat0, ONLY : init_etat0, etat0
[4855]	43	USE icolmdz_param_gravity_wave, ONLY: init_param_gravity_wave
	44	IMPLICIT NONE
	45	!$OMP PARALLEL
	46	init_etat0_plugin => init_etat0
	47	etat0_plugin => etat0
	48	CALL init_param_gravity_wave
	49	!$OMP END PARALLEL
	50	END SUBROUTINE pre_initialize_physics
	51
[4498]	52	SUBROUTINE initialize_physics
	53	USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
	54	! from dynamico
	55	USE domain_mod
	56	USE dimensions
	57	USE mpi_mod
	58	USE mpipara
	59	USE disvert_mod
	60	USE xios_mod
	61	USE time_mod , init_time_icosa=> init_time
	62	USE transfert_mod
	63	USE nudging_mod, ONLY : lam_halo_scheme
	64
	65	! from LMDZ
	66	USE mod_grid_phy_lmdz, ONLY : unstructured
	67	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
	68	USE transfert_mod
	69	USE physics_distribution_mod, ONLY : init_physics_distribution
[4501]	70	USE infotrac_phy, ONLY: init_infotrac_phy
[4855]	71	USE icolmdz_param_gravity_wave, ONLY: param_gravity_wave
[4501]	72
[4498]	73
	74
	75	IMPLICIT NONE
	76	INTEGER :: ind,i,j,ij,pos,h
	77	REAL(rstd),POINTER :: bounds_lon(:,:)
	78	REAL(rstd),POINTER :: bounds_lat(:,:)
	79
	80	REAL(rstd),ALLOCATABLE :: latfi(:)
	81	REAL(rstd),ALLOCATABLE :: lonfi(:)
	82	REAL(rstd),ALLOCATABLE :: airefi(:)
	83	REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
	84	REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
	85	LOGICAL ,ALLOCATABLE :: outside(:,:)
	86	LOGICAL ,ALLOCATABLE :: outside_tmp(:,:)
	87	LOGICAL ,POINTER :: out(:,:)
	88	! REAL(rstd) :: pseudoalt(llm)
	89
	90	INTEGER :: nbp_phys, nbp_phys_glo
	91
	92	!$OMP PARALLEL
	93	CALL allocate_field(f_bounds_lon,field_t,type_real,6)
	94	CALL allocate_field(f_bounds_lat,field_t,type_real,6)
	95	CALL allocate_field(f_p,field_t,type_real,llm+1,name="p_in")
	96	CALL allocate_field(f_pks,field_t,type_real)
	97	CALL allocate_field(f_pk,field_t,type_real,llm)
	98	CALL allocate_field(f_p_layer,field_t,type_real,llm,name="p_layer_in")
	99	CALL allocate_field(f_theta,field_t,type_real,llm)
	100	CALL allocate_field(f_phi,field_t,type_real,llm,name="phi_in")
	101	CALL allocate_field(f_Temp,field_t,type_real,llm,name="Temp_in")
	102	CALL allocate_field(f_ulon,field_t,type_real,llm,name="ulon_in")
	103	CALL allocate_field(f_ulat,field_t,type_real,llm,name="ulat_in")
	104	CALL allocate_field(f_vort,field_z,type_real,llm,name="vort_in")
	105	CALL allocate_field(f_vortc,field_t,type_real,llm,name="vortc_in")
	106	CALL allocate_field(f_dulon,field_t,type_real,llm,name="dulon_out")
	107	CALL allocate_field(f_dulat,field_t,type_real,llm,name="dulat_out")
	108	CALL allocate_field(f_dTemp,field_t,type_real,llm,name="dTemp_out")
	109	CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,name="dq_out")
	110	CALL allocate_field(f_dps,field_t,type_real,name="dps_out")
	111	CALL allocate_field(f_duc,field_t,type_real,3,llm)
	112
	113	CALL init_message(f_dps,req_i0,req_dps0)
	114	CALL init_message(f_dulon,req_i0,req_dulon0)
	115	CALL init_message(f_dulat,req_i0,req_dulat0)
	116	CALL init_message(f_dTemp,req_i0,req_dTemp0)
	117	CALL init_message(f_dq,req_i0,req_dq0)
	118	!$OMP END PARALLEL
	119
	120	nbp_phys=0
	121	DO ind=1,ndomain
	122	CALL swap_dimensions(ind)
	123
	124	ALLOCATE(outside(ii_begin:ii_end,jj_begin:jj_end)) ! for limited area : don't take cells arround the border
	125	ALLOCATE(outside_tmp(ii_begin-1:ii_end+1,jj_begin-1:jj_end+1)) ! for limited area : don't take cells arround the border
	126	out=>domain(ind)%outside
	127	DO j=jj_begin,jj_end
	128	DO i=ii_begin,ii_end
	129	outside(i,j)= out(i+1,j) .OR. & ! right
	130	out(i,j+1 ) .OR. & ! rup
	131	out(i-1 ,j+1) .OR. & ! lup
	132	out(i-1 ,j) .OR. & !left
	133	out(i ,j-1) .OR. & !ldown
	134	out(i+1,j-1) !rdown
	135	ENDDO
	136	ENDDO
	137
	138	outside_tmp=.FALSE.
	139	outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside
	140
	141	DO h=1,lam_halo_scheme-1 ! do not compute physic on limited area halo
	142	DO j=jj_begin,jj_end
	143	DO i=ii_begin,ii_end
	144	outside(i,j) = outside_tmp(i,j) .OR. &
	145	outside_tmp(i+1,j) .OR. & ! right
	146	outside_tmp(i,j+1 ) .OR. & ! rup
	147	outside_tmp(i-1 ,j+1) .OR. & ! lup
	148	outside_tmp(i-1 ,j) .OR. & !left
	149	outside_tmp(i ,j-1) .OR. & !ldown
	150	outside_tmp(i+1,j-1) !rdown
	151	ENDDO
	152	ENDDO
	153	outside_tmp(ii_begin:ii_end,jj_begin:jj_end)=outside
	154	ENDDO
	155
	156	DO j=jj_begin,jj_end
	157	DO i=ii_begin,ii_end
	158	IF (domain(ind)%own(i,j) .AND. .NOT.outside(i,j)) nbp_phys=nbp_phys+1
	159	ENDDO
	160	ENDDO
	161	DEALLOCATE(outside)
	162	DEALLOCATE(outside_tmp)
	163	ENDDO
	164
	165
	166	!initialize LMDZ5 physic mpi decomposition
	167	CALL MPI_ALLREDUCE(nbp_phys,nbp_phys_glo,1,MPI_INTEGER,MPI_SUM,comm_icosa,ierr)
	168	CALL init_physics_distribution(unstructured, 6, nbp_phys, 1, nbp_phys_glo, llm, comm_icosa)
	169
	170	DO ind=1,ndomain
	171	CALL swap_dimensions(ind)
	172	CALL swap_geometry(ind)
	173	bounds_lon=f_bounds_lon(ind)
	174	bounds_lat=f_bounds_lat(ind)
	175	DO j=jj_begin,jj_end
	176	DO i=ii_begin,ii_end
	177	ij=(j-1)*iim+i
	178	CALL xyz2lonlat(xyz_v(ij+z_rup,:), bounds_lon(ij,1), bounds_lat(ij,1))
	179	CALL xyz2lonlat(xyz_v(ij+z_up,:), bounds_lon(ij,2), bounds_lat(ij,2))
	180	CALL xyz2lonlat(xyz_v(ij+z_lup,:), bounds_lon(ij,3), bounds_lat(ij,3))
	181	CALL xyz2lonlat(xyz_v(ij+z_ldown,:), bounds_lon(ij,4), bounds_lat(ij,4))
	182	CALL xyz2lonlat(xyz_v(ij+z_down,:), bounds_lon(ij,5), bounds_lat(ij,5))
	183	CALL xyz2lonlat(xyz_v(ij+z_rdown,:), bounds_lon(ij,6), bounds_lat(ij,6))
	184	ENDDO
	185	ENDDO
	186	ENDDO
[4501]	187
	188	CALL init_infotrac_phy
	189
[4498]	190
	191	!$OMP PARALLEL
	192	CALL initialize_physics_omp
[4855]	193	CALL param_gravity_wave
[4498]	194	!$OMP END PARALLEL
	195
	196	CALL xios_set_context
	197
	198	END SUBROUTINE initialize_physics
	199
	200
	201	SUBROUTINE initialize_physics_omp
	202	USE distrib_icosa_lmdz_mod, ONLY : init_distrib_icosa_lmdz, transfer_icosa_to_lmdz
	203	! from dynamico
	204	USE domain_mod
	205	USE dimensions
	206	USE mpi_mod
	207	USE mpipara
	208	USE disvert_mod
	209	USE earth_const, ONLY: scale_height
	210	USE xios_mod
	211	USE time_mod , init_time_icosa=> init_time
	212	USE omp_para
	213
	214	! from LMDZ
	215	USE mod_grid_phy_lmdz, ONLY : unstructured, klon_glo
	216	USE mod_phys_lmdz_para, ONLY: klon_omp, reduce_min_lmdz => reduce_min , gather_lmdz => gather , bcast_lmdz => bcast
	217	USE time_phylmdz_mod, ONLY: init_time_lmdz => init_time
	218	USE transfert_mod
	219	! USE physics_distribution_mod, ONLY : init_physics_distribution
	220	USE geometry_mod, ONLY : init_geometry
	221	USE vertical_layers_mod, ONLY : init_vertical_layers
	222	USE infotrac_phy, ONLY : init_infotrac_phy
	223	USE inifis_mod, ONLY : inifis
[4501]	224	USE readTracFiles_mod, ONLY: trac_type, isot_type
[4498]	225	USE tracer_icosa_mod, ONLY : tracs
	226	USE readTracFiles_mod, ONLY: delPhase
	227	! USE phyaqua_mod, ONLY : iniaqua
[5018]	228	USE phystokenc_mod, ONLY : init_phystokenc ! use to generate mass flow files for offline mode
[4498]	229
	230	IMPLICIT NONE
	231
	232
	233
	234	INTEGER :: ind,i,j,k,ij,pos
	235	REAL(rstd),POINTER :: bounds_lon(:,:)
	236	REAL(rstd),POINTER :: bounds_lat(:,:)
	237
	238	REAL(rstd),ALLOCATABLE :: latfi(:)
	239	REAL(rstd),ALLOCATABLE :: lonfi(:)
	240	REAL(rstd),ALLOCATABLE :: airefi(:)
	241	REAL(rstd),ALLOCATABLE :: bounds_latfi(:,:)
	242	REAL(rstd),ALLOCATABLE :: bounds_lonfi(:,:)
	243	REAL(rstd),ALLOCATABLE :: ind_cell_glo_r(:)
	244	INTEGER, ALLOCATABLE :: ind_cell_glo(:)
	245	INTEGER, ALLOCATABLE :: ind_cell_glo_tot(:)
	246	INTEGER, ALLOCATABLE :: cell_glo_tot(:)
	247	INTEGER :: ncell_glo_tot
	248
	249	REAL(rstd) :: pseudoalt(llm)
	250	REAL(rstd) :: aps(llm)
	251	REAL(rstd) :: bps(llm)
	252	REAL(rstd) :: scaleheight
	253
	254	INTEGER :: run_length
	255	REAL :: day_length ! length of a day (s) ! SAVEd to be OpenMP shared <--- NO!!!!
	256	INTEGER :: annee_ref
	257	INTEGER :: day_ref
	258	INTEGER :: day_ini
	259	REAL :: start_time
	260	REAL :: physics_timestep
	261
	262	! Tracer stuff (SAVEd when needed to be OpenMP shared)
	263	INTEGER :: nq
	264	INTEGER :: nqo, nbtr, nbtr_inca
	265	CHARACTER(len=256) :: type_trac
	266	INTEGER,ALLOCATABLE :: conv_flg(:) ! conv_flg(it)=0 : convection desactivated for tracer number it
	267	INTEGER,ALLOCATABLE :: pbl_flg(:) ! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it
	268	CHARACTER(len=8),ALLOCATABLE :: solsym(:) ! tracer name from inca
	269
	270	TYPE(t_field),POINTER,SAVE :: f_ind_cell_glo(:)
	271
	272	INTEGER :: iflag_phys
	273
	274	INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca ! index of horizontal trasport schema
	275	INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca ! index of vertical trasport schema
	276
	277	TYPE(trac_type) :: tracers_ico2lmdz(nqtot) !=== TRACERS DESCRIPTORS VECTOR
	278	TYPE(isot_type) :: isotopes_ico2lmdz(1) !=== ISOTOPES PARAMETERS VECTOR
	279	INTEGER :: iq
	280
	281	CHARACTER(LEN=3) :: descrq(30) !--- Advection scheme description tags
	282	logical, save :: first = .TRUE.
	283
[5018]	284	LOGICAL :: offline
	285	INTEGER :: offline_time
[4498]	286
[5018]	287
[4498]	288	CALL init_distrib_icosa_lmdz
	289
	290	ALLOCATE(latfi(klon_omp))
	291	ALLOCATE(lonfi(klon_omp))
	292	ALLOCATE(airefi(klon_omp))
	293	ALLOCATE(bounds_latfi(klon_omp,6))
	294	ALLOCATE(bounds_lonfi(klon_omp,6))
	295	ALLOCATE(ind_cell_glo_r(klon_omp))
	296	ALLOCATE(ind_cell_glo(klon_omp))
	297
	298	CALL transfer_icosa_to_lmdz(geom%lat_i,latfi)
	299	CALL transfer_icosa_to_lmdz(geom%lon_i,lonfi)
	300	CALL transfer_icosa_to_lmdz(f_bounds_lat,bounds_latfi)
	301	CALL transfer_icosa_to_lmdz(f_bounds_lon,bounds_lonfi)
	302	CALL transfer_icosa_to_lmdz(geom%Ai,airefi)
	303
	304	CALL allocate_field(f_ind_cell_glo,field_t,type_real)
	305
	306	DO ind=1,ndomain
	307	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
	308	CALL swap_dimensions(ind)
	309	CALL swap_geometry(ind)
	310	DO j=jj_begin,jj_end
	311	DO i=ii_begin,ii_end
	312	ij=(j-1)*iim+i
	313	f_ind_cell_glo(ind)%rval2d(ij)=domain(ind)%assign_cell_glo(i,j)
	314	ENDDO
	315	ENDDO
	316	ENDDO
	317
	318
	319	CALL transfer_icosa_to_lmdz(f_ind_cell_glo,ind_cell_glo_r)
	320	CALL deallocate_field(f_ind_cell_glo)
	321	ind_cell_glo=INT(ind_cell_glo_r)
	322	DEALLOCATE(ind_cell_glo_r)
	323
	324
	325	CALL reduce_min_lmdz(MINVAL(-ind_cell_glo),ncell_glo_tot) ! reduce_max does not exist in lmdz, use reduce_min
	326	CALL bcast_lmdz(ncell_glo_tot)
	327	ncell_glo_tot=-ncell_glo_tot
	328	ALLOCATE(cell_glo_tot(0:ncell_glo_tot))
	329	ALLOCATE(ind_cell_glo_tot(klon_glo))
	330	CALL gather_lmdz(ind_cell_glo,ind_cell_glo_tot)
	331	CALL bcast_lmdz(ind_cell_glo_tot)
	332
	333	cell_glo_tot=-1
	334	DO i=1,klon_glo
	335	cell_glo_tot(ind_cell_glo_tot(i))= 0
	336	ENDDO
	337
	338	pos=0
	339	DO i=0,ncell_glo_tot
	340	IF (cell_glo_tot(i)/=-1) THEN
	341	cell_glo_tot(i) = pos
	342	pos=pos + 1
	343	ENDIF
	344	ENDDO
	345
	346	DO i=1,klon_omp
	347	ind_cell_glo(i)=cell_glo_tot(ind_cell_glo(i))
	348	ENDDO
	349
	350	ind_cell_glo = ind_cell_glo + 1 ! lmdz expect global indices begining to 1 not 0
	351
	352
	353	! CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_mpi,1,MPI_INTEGER,comm_icosa,ierr)
	354	!
	355	! displ(0)=0
	356	! DO i=1,mpi_size-1
	357	! displ(i)=displ(i-1)+ncell_mpi(i-1)
	358	! ENDDO
	359	!
	360	! ALLOCATE(ind_glo_tot(ncell_tot))
	361	! ALLOCATE(cell_glo_tot(0:ncell_glo_tot-1))
	362	!
	363	! cell_glo_tot(:)= -1
	364	! CALL MPI_ALLGATHERV(ind_glo, ncell, MPI_INTEGER, ind_glo_tot, ncell_mpi, displ, MPI_INTEGER, comm_icosa,ierr)
	365	!
	366	! DO i=1,ncell_tot
	367	! cell_glo_tot(ind_glo_tot(i))= 0
	368	! ENDDO
	369	!
	370	! ncell_glo=0
	371	! DO i=0,ncell_glo_tot-1
	372	! IF (cell_glo_tot(i)/=-1) THEN
	373	! cell_glo_tot(i) = ncell_glo
	374	! ncell_glo=ncell_glo + 1
	375	! ENDIF
	376	! ENDDO
	377	!
	378	! DO i=1,ncell
	379	! ind_glo(i)=cell_glo_tot(ind_glo(i))
	380	! ENDDO
	381
	382
	383	CALL init_geometry(klon_omp,lonfi, latfi, bounds_lonfi, bounds_latfi, airefi, ind_cell_glo)
	384
	385	scaleheight=scale_height/1000. ! Atmospheric scale height (km)
	386	aps(1:llm)=0.5*(ap(1:llm)+ap(2:llm+1))
	387	bps(1:llm)=0.5*(bp(1:llm)+bp(2:llm+1))
	388	pseudoalt(:)=-scaleheight*log(presnivs(:)/preff)
	389	CALL init_vertical_layers(llm,preff,scaleheight,ap,bp,aps,bps,presnivs,presinter,pseudoalt)
	390
	391	! Initialize tracer names, numbers, etc. for physics
	392	!Config Key = type_trac
	393	!Config Desc = Choix de couplage avec model de chimie INCA ou REPROBUS
	394	!Config Def = lmdz
	395	!Config Help =
	396	!Config 'lmdz' = pas de couplage, pur LMDZ
	397	!Config 'lmdz\|inca' = model de chime INCA
	398	!Config 'lmdz\|repr' = model de chime REPROBUS
	399	type_trac = 'lmdz'
	400	CALL getin('type_trac',type_trac)
	401
	402	descrq( 1: 2) = ['LMV','BAK']
	403	descrq(10:20) = ['VL1','VLP','FH1','FH2','VLH',' ','PPM','PPS','PPP',' ','SLP']
	404	descrq(30) = 'PRA'
	405
	406	nqo = 0
	407	DO iq=1,nqtot
	408
	409	tracers_ico2lmdz(iq)%name = tracs(iq)%name
	410
	411	tracers_ico2lmdz(iq)%gen0Name = tracs(iq)%name
	412	tracers_ico2lmdz(iq)%phase = tracs(iq)%phase
	413	tracers_ico2lmdz(iq)%iadv = tracs(iq)%iadv
	414
	415	IF (tracs(iq)%component .eq. "dynamico") then
	416	tracers_ico2lmdz(iq)%component='lmdz'
	417	ELSE
	418	tracers_ico2lmdz(iq)%component=tracs(iq)%component
	419	ENDIF
	420
	421	if (tracers_ico2lmdz(iq)%iadv .ne. 0 ) tracers_ico2lmdz(iq)%isAdvected=.true.
	422
	423	tracers_ico2lmdz(iq)%longName = tracers_ico2lmdz(iq)%name
	424	IF(tracers_ico2lmdz(iq)%iadv > 0) tracers_ico2lmdz(iq)%longName=TRIM(tracers_ico2lmdz(iq)%name)//descrq(tracers_ico2lmdz(iq)%iadv)
	425
	426	tracers_ico2lmdz(iq)%isInPhysics= delPhase(tracers_ico2lmdz(iq)%gen0Name) /= 'H2O' .OR. tracers_ico2lmdz(iq)%component /= 'lmdz'
	427	tracers_ico2lmdz(iq)%iGeneration = 0
	428
	429	ENDDO
	430
	431	nqo = COUNT(delPhase(tracs(:)%name) == 'H2O' .AND. tracers_ico2lmdz(:)%component == 'lmdz') !--- Number of water phases
	432
	433	isotopes_ico2lmdz(1)%parent='H2O'
	434	isotopes_ico2lmdz(1)%phase='gls'
	435	isotopes_ico2lmdz(1)%nphas=3
	436
	437
	438	nbtr=nqtot-nqo
	439
	440	ALLOCATE(conv_flg(nbtr))
	441	ALLOCATE(pbl_flg(nbtr))
	442
	443	conv_flg(:) = 1 ! convection activated for all tracers
	444	pbl_flg(:) = 1 ! boundary layer activated for all tracers
	445
	446
[4501]	447	CALL init_infotrac_phy
[4498]	448
	449	! Initialize physical constant
	450	day_length=86400
	451	CALL getin('day_length',day_length)
	452	CALL inifis(day_length,radius,g,kappa*cpp,cpp)
	453
	454
	455
	456	! init time
	457	annee_ref=2015
	458	CALL getin("anneeref",annee_ref)
	459
	460	day_ref=1
	461	CALL getin("dayref",day_ref)
	462
	463	physics_timestep=dt*itau_physics
	464	run_length=itaumax*dt
	465	ndays=NINT(run_length/day_length)
	466
	467	day_ini=INT(itau0*dt/day_length)+day_ref
	468	start_time= itau0dt/day_length-INT(itau0dt/day_length)
	469
	470	CALL init_time_lmdz(annee_ref, day_ref, day_ini, start_time, int(ndays), physics_timestep)
	471
[5018]	472
	473	! Init Offline mode
	474	offline = .FALSE.
	475	CALL getin('offline',offline)
	476
	477	! Choosing storage frequencies for offline mass flow files
	478	! offline_time=12 2h=1day/12
	479	! offline_time=8 3h=1day/8
	480	offline_time = 8
	481	CALL getin('offline_time',offline_time)
	482
	483	! Copy over "offline" settings
	484	! Flag and number of time steps for flux calculation and output
	485	CALL init_phystokenc(offline,int(day_length/(offline_time*physics_timestep)))
	486
	487
[4498]	488	! Additional initializations for aquaplanets
	489	! CALL getin("iflag_phys",iflag_phys)
	490	! IF (iflag_phys>=100) THEN
	491	! CALL iniaqua(klon_omp, iflag_phys)
	492	! END IF
	493
	494
	495
	496	#ifdef INCA
	497	CONTAINS
	498
	499	SUBROUTINE init_chem_trac()
	500	IMPLICIT NONE
	501
	502	CALL Init_chem_inca_trac(nbtr)
	503
	504	END SUBROUTINE init_chem_trac
	505
	506	SUBROUTINE init_chem_transport()
	507
	508	IMPLICIT NONE
	509
	510	CALL init_transport(solsym, conv_flg,pbl_flg, hadv_inca, vadv_inca)
	511
	512	END SUBROUTINE init_chem_transport
	513
	514
	515	#else
	516	CONTAINS
	517	SUBROUTINE init_chem_trac()
	518	IMPLICIT NONE
	519
	520	END SUBROUTINE init_chem_trac
	521
	522	SUBROUTINE init_chem_transport()
	523
	524	IMPLICIT NONE
	525
	526	END SUBROUTINE init_chem_transport
	527
	528
	529	#endif
	530
	531
	532
	533	END SUBROUTINE initialize_physics_omp
	534
	535
	536
	537
	538	SUBROUTINE physics
	539	USE icosa
	540	USE time_mod
	541	USE disvert_mod
	542	USE transfert_mod
	543	USE mpipara
	544	USE xios_mod
	545	USE wxios
	546	USE trace
	547	USE distrib_icosa_lmdz_mod, ONLY : transfer_icosa_to_lmdz, transfer_lmdz_to_icosa
	548	USE physics_external_mod, ONLY : it, f_phis, f_ps, f_theta_rhodz, f_u, f_wflux, f_q
	549	USE write_field_mod
	550	USE checksum_mod
	551	USE vorticity_mod
	552
	553	! from LMDZ
	554	USE mod_phys_lmdz_omp_data, ONLY: klon_omp
	555	USE geometry_mod, ONLY : cell_area
	556	USE physiq_mod, ONLY: physiq
[4855]	557	USE icolmdz_param_gravity_wave, ONLY: param_gravity_wave
[4498]	558	IMPLICIT NONE
	559
	560	REAL(rstd),POINTER :: phis(:)
	561	REAL(rstd),POINTER :: ps(:)
	562	REAL(rstd),POINTER :: theta_rhodz(:,:,:)
	563	REAL(rstd),POINTER :: u(:,:)
	564	REAL(rstd),POINTER :: wflux(:,:)
	565	REAL(rstd),POINTER :: q(:,:,:)
	566	REAL(rstd),POINTER :: p(:,:)
	567	REAL(rstd),POINTER :: pks(:)
	568	REAL(rstd),POINTER :: pk(:,:)
	569	REAL(rstd),POINTER :: p_layer(:,:)
	570	REAL(rstd),POINTER :: theta(:,:)
	571	REAL(rstd),POINTER :: phi(:,:)
	572	REAL(rstd),POINTER :: Temp(:,:)
	573	REAL(rstd),POINTER :: ulon(:,:)
	574	REAL(rstd),POINTER :: ulat(:,:)
	575	REAL(rstd),POINTER :: vort(:,:)
	576	REAL(rstd),POINTER :: vortc(:,:)
	577	REAL(rstd),POINTER :: dulon(:,:)
	578	REAL(rstd),POINTER :: dulat(:,:)
	579	REAL(rstd),POINTER :: dTemp(:,:)
	580	REAL(rstd),POINTER :: dq(:,:,:)
	581	REAL(rstd),POINTER :: dps(:)
	582	REAL(rstd),POINTER :: duc(:,:,:)
	583
	584
	585	INTEGER :: ind,l
	586
	587	REAL(rstd),ALLOCATABLE,SAVE :: ps_phy(:)
	588	!$OMP THREADPRIVATE(ps_phy)
	589	REAL(rstd),ALLOCATABLE,SAVE :: p_phy(:,:)
	590	!$OMP THREADPRIVATE(p_phy)
	591	REAL(rstd),ALLOCATABLE,SAVE :: p_layer_phy(:,:)
	592	!$OMP THREADPRIVATE(p_layer_phy)
	593	REAL(rstd),ALLOCATABLE,SAVE :: Temp_phy(:,:)
	594	!$OMP THREADPRIVATE(Temp_phy)
	595	REAL(rstd),ALLOCATABLE,SAVE :: phis_phy(:)
	596	!$OMP THREADPRIVATE(phis_phy)
	597	REAL(rstd),ALLOCATABLE,SAVE :: phi_phy(:,:)
	598	!$OMP THREADPRIVATE(phi_phy)
	599	REAL(rstd),ALLOCATABLE,SAVE :: ulon_phy(:,:)
	600	!$OMP THREADPRIVATE(ulon_phy)
	601	REAL(rstd),ALLOCATABLE,SAVE :: ulat_phy(:,:)
	602	!$OMP THREADPRIVATE(ulat_phy)
	603	REAL(rstd),ALLOCATABLE,SAVE :: rot_phy(:,:)
	604	!$OMP THREADPRIVATE(rot_phy)
	605	REAL(rstd),ALLOCATABLE,SAVE :: q_phy(:,:,:)
	606	!$OMP THREADPRIVATE(q_phy)
	607	REAL(rstd),ALLOCATABLE,SAVE :: wflux_phy(:,:)
	608	!$OMP THREADPRIVATE(wflux_phy)
	609	REAL(rstd),ALLOCATABLE,SAVE :: dulon_phy(:,:)
	610	!$OMP THREADPRIVATE(dulon_phy)
	611	REAL(rstd),ALLOCATABLE,SAVE :: dulat_phy(:,:)
	612	!$OMP THREADPRIVATE(dulat_phy)
	613	REAL(rstd),ALLOCATABLE,SAVE :: dTemp_phy(:,:)
	614	!$OMP THREADPRIVATE(dTemp_phy)
	615	REAL(rstd),ALLOCATABLE,SAVE :: dq_phy(:,:,:)
	616	!$OMP THREADPRIVATE(dq_phy)
	617	REAL(rstd),ALLOCATABLE,SAVE :: dps_phy(:)
	618	!$OMP THREADPRIVATE(dps_phy)
	619	REAL(rstd) :: dtphy
	620	LOGICAL :: debut
	621	LOGICAL :: lafin
	622	LOGICAL,SAVE :: first=.TRUE.
	623	!$OMP THREADPRIVATE(first)
	624
	625
	626	IF(first) THEN
	627	debut=.TRUE.
	628	ELSE
	629	debut=.FALSE.
	630	ENDIF
	631
	632
	633	IF(it-itau0>=itaumax) THEN
	634	lafin=.TRUE.
	635	ELSE
	636	lafin=.FALSE.
	637	ENDIF
	638
	639	IF (first) THEN
	640	first=.FALSE.
	641	CALL init_message(f_u,req_e1_vect,req_u)
	642	CALL init_message(f_vort,req_z1_scal,req_z)
	643	ALLOCATE(ps_phy(klon_omp))
	644	ALLOCATE(p_phy(klon_omp,llm+1))
	645	ALLOCATE(p_layer_phy(klon_omp,llm))
	646	ALLOCATE(Temp_phy(klon_omp,llm))
	647	ALLOCATE(phis_phy(klon_omp))
	648	ALLOCATE(phi_phy(klon_omp,llm))
	649	ALLOCATE(ulon_phy(klon_omp,llm))
	650	ALLOCATE(ulat_phy(klon_omp,llm))
	651	ALLOCATE(rot_phy(klon_omp,llm))
	652	ALLOCATE(q_phy(klon_omp,llm,nqtot))
	653	ALLOCATE(wflux_phy(klon_omp,llm))
	654	ALLOCATE(dulon_phy(klon_omp,llm))
	655	ALLOCATE(dulat_phy(klon_omp,llm))
	656	ALLOCATE(dTemp_phy(klon_omp,llm))
	657	ALLOCATE(dq_phy(klon_omp,llm,nqtot))
	658	ALLOCATE(dps_phy(klon_omp))
	659	!$OMP BARRIER
[4855]	660
	661	! CALL param_gravity_wave
[4498]	662	ENDIF
	663
	664
	665	!$OMP MASTER
	666	! CALL update_calendar(it)
	667	!$OMP END MASTER
	668	!$OMP BARRIER
	669	dtphy=itau_physics*dt
	670
	671
	672
	673	CALL transfert_message(f_u,req_u)
	674	DO ind=1,ndomain
	675	IF (assigned_domain(ind)) THEN
	676	CALL swap_dimensions(ind)
	677	CALL swap_geometry(ind)
	678	u=f_u(ind)
	679	vort=f_vort(ind)
	680	CALL compute_vorticity(u,vort)
	681	ENDIF
	682	ENDDO
	683
	684	CALL transfert_message(f_vort,req_z)
	685
	686
	687	DO ind=1,ndomain
	688	CALL swap_dimensions(ind)
	689	IF (assigned_domain(ind)) THEN
	690	CALL swap_geometry(ind)
	691
	692	phis=f_phis(ind)
	693	ps=f_ps(ind)
	694	theta_rhodz=f_theta_rhodz(ind)
	695	u=f_u(ind)
	696	q=f_q(ind)
	697	wflux=f_wflux(ind)
	698	p=f_p(ind)
	699	pks=f_pks(ind)
	700	pk=f_pk(ind)
	701	p_layer=f_p_layer(ind)
	702	theta=f_theta(ind)
	703	phi=f_phi(ind)
	704	Temp=f_Temp(ind)
	705	ulon=f_ulon(ind)
	706	ulat=f_ulat(ind)
	707	vort=f_vort(ind)
	708	vortc=f_vortc(ind)
	709
	710	CALL grid_icosa_to_physics
	711
	712	ENDIF
	713	ENDDO
	714
	715	!$OMP BARRIER
	716	!$OMP MASTER
	717	CALL SYSTEM_CLOCK(start_clock)
	718	!$OMP END MASTER
	719	CALL trace_start("physic")
	720	! CALL trace_off()
	721
	722
	723	! CALL writeField("p_in",f_p)
	724	! CALL writeField("p_layer_in",f_p_layer)
	725	! CALL writeField("phi_in",f_phi)
	726	! CALL writeField("phis_in",f_phis)
	727	! CALL writeField("ulon_in",f_ulon)
	728	! CALL writeField("ulat_in",f_ulat)
	729	! CALL writeField("Temp_in",f_Temp)
	730	! CALL writeField("q_in",f_q)
	731	! CALL writeField("wflux_in",f_wflux)
	732	! CALL writeField("vortc",f_vortc)
	733
	734	! CALL checksum(f_p)
	735	! CALL checksum(f_p_layer)
	736	! CALL checksum(f_phi)
	737	! CALL checksum(f_phis)
	738	! CALL checksum(f_ulon)
	739	! CALL checksum(f_ulat)
	740	! CALL checksum(f_Temp)
	741	! CALL checksum(f_q)
	742	! CALL checksum(f_wflux)
	743
	744	CALL transfer_icosa_to_lmdz(f_p , p_phy)
	745	CALL transfer_icosa_to_lmdz(f_p_layer, p_layer_phy)
	746	CALL transfer_icosa_to_lmdz(f_phi , phi_phy)
	747	CALL transfer_icosa_to_lmdz(f_phis , phis_phy )
	748	CALL transfer_icosa_to_lmdz(f_ulon , ulon_phy )
	749	CALL transfer_icosa_to_lmdz(f_ulat , ulat_phy)
	750	CALL transfer_icosa_to_lmdz(f_vortc , rot_phy)
	751	CALL transfer_icosa_to_lmdz(f_Temp , Temp_phy)
	752	CALL transfer_icosa_to_lmdz(f_q , q_phy)
	753	CALL transfer_icosa_to_lmdz(f_wflux , wflux_phy)
	754
	755	DO l=1,llm
	756	wflux_phy(:,l) = - wflux_phy(:,l)*cell_area(:)
	757	phi_phy(:,l)=phi_phy(:,l)-phis_phy(:)
	758	ENDDO
	759
	760	CALL wxios_set_context()
	761
	762	! Ehouarn: rot_phy() not implemented!! Set it to zero for now
	763	! rot_phy(:,:)=0
	764	CALL physiq(klon_omp, llm, debut, lafin, dtphy, &
	765	p_phy, p_layer_phy, phi_phy, phis_phy, presnivs, &
	766	ulon_phy, ulat_phy, rot_phy, Temp_phy, q_phy, wflux_phy, &
	767	dulon_phy, dulat_phy, dTemp_phy, dq_phy, dps_phy)
	768
	769	CALL transfer_lmdz_to_icosa(dulon_phy, f_dulon )
	770	CALL transfer_lmdz_to_icosa(dulat_phy, f_dulat )
	771	CALL transfer_lmdz_to_icosa(dTemp_phy, f_dTemp )
	772	CALL transfer_lmdz_to_icosa(dq_phy , f_dq )
	773	CALL transfer_lmdz_to_icosa(dps_phy , f_dps )
	774
	775	! CALL writeField("dulon_out",f_dulon)
	776	! CALL writeField("dulat_out",f_dulat)
	777	! CALL writeField("dTemp_out",f_dTemp)
	778	! CALL writeField("dq_out",f_dq)
	779	! CALL writeField("dps_out",f_dps)
	780
	781	! CALL checksum(f_dulon)
	782	! CALL checksum(f_dulat)
	783	! CALL checksum(f_dTemp)
	784	! CALL checksum(f_dq)
	785	! CALL checksum(f_dps)
	786
	787	CALL send_message(f_dps,req_dps0)
	788	CALL send_message(f_dulon,req_dulon0)
	789	CALL send_message(f_dulat,req_dulat0)
	790	CALL send_message(f_dTemp,req_dTemp0)
	791	CALL send_message(f_dq,req_dq0)
	792
	793	CALL wait_message(req_dps0)
	794	CALL wait_message(req_dulon0)
	795	CALL wait_message(req_dulat0)
	796	CALL wait_message(req_dTemp0)
	797	CALL wait_message(req_dq0)
	798
	799
	800	! CALL trace_on()
	801	CALL trace_end("physic")
	802	!$OMP MASTER
	803	CALL SYSTEM_CLOCK(stop_clock)
	804	count_clock=count_clock+stop_clock-start_clock
	805	!$OMP END MASTER
	806
	807	!$OMP BARRIER
	808
	809	DO ind=1,ndomain
	810	CALL swap_dimensions(ind)
	811	IF (assigned_domain(ind)) THEN
	812	CALL swap_geometry(ind)
	813
	814	theta_rhodz=f_theta_rhodz(ind)
	815	u=f_u(ind)
	816	q=f_q(ind)
	817	ps=f_ps(ind)
	818	dulon=f_dulon(ind)
	819	dulat=f_dulat(ind)
	820	Temp=f_temp(ind)
	821	dTemp=f_dTemp(ind)
	822	dq=f_dq(ind)
	823	dps=f_dps(ind)
	824	duc=f_duc(ind)
	825	p=f_p(ind)
	826	pks=f_pks(ind)
	827	pk=f_pk(ind)
	828
	829	CALL grid_physics_to_icosa
	830	ENDIF
	831	ENDDO
	832
	833	!$OMP BARRIER
	834	CALL xios_set_context
	835
	836
	837	CONTAINS
	838
	839	SUBROUTINE grid_icosa_to_physics
	840	USE pression_mod
	841	USE exner_mod
	842	USE theta2theta_rhodz_mod
	843	USE geopotential_mod
[4857]	844	USE wind_from_lonlat_mod
[4498]	845	USE omp_para
	846	IMPLICIT NONE
	847
	848	REAL(rstd) :: uc(3)
	849	INTEGER :: i,j,ij,l
	850
	851	! compute pression
	852
	853	DO l = ll_begin,ll_endp1
	854	DO j=jj_begin,jj_end
	855	DO i=ii_begin,ii_end
	856	ij=(j-1)*iim+i
	857	p(ij,l) = ap(l) + bp(l) * ps(ij)
	858	ENDDO
	859	ENDDO
	860	ENDDO
	861
	862	!$OMP BARRIER
	863
	864	! compute exner
	865
	866	IF (is_omp_first_level) THEN
	867	DO j=jj_begin,jj_end
	868	DO i=ii_begin,ii_end
	869	ij=(j-1)*iim+i
	870	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
	871	ENDDO
	872	ENDDO
	873	ENDIF
	874
	875	! 3D : pk
	876	DO l = ll_begin,ll_end
	877	DO j=jj_begin,jj_end
	878	DO i=ii_begin,ii_end
	879	ij=(j-1)*iim+i
	880	pk(ij,l) = cpp * ((.5/preff)(p(ij,l)+p(ij,l+1))) * kappa
	881	ENDDO
	882	ENDDO
	883	ENDDO
	884
	885	!$OMP BARRIER
	886
	887	! compute theta, temperature and pression at layer
	888	DO l = ll_begin, ll_end
	889	DO j=jj_begin,jj_end
	890	DO i=ii_begin,ii_end
	891	ij=(j-1)*iim+i
	892	theta(ij,l) = theta_rhodz(ij,l,1) / ((p(ij,l)-p(ij,l+1))/g)
	893	Temp(ij,l) = theta(ij,l) * pk(ij,l) / cpp
	894	p_layer(ij,l)=preff(pk(ij,l)/cpp)*(1./kappa)
	895	ENDDO
	896	ENDDO
	897	ENDDO
	898
	899
	900	!!! Compute geopotential
	901
	902	! for first layer
	903	IF (is_omp_first_level) THEN
	904	DO j=jj_begin,jj_end
	905	DO i=ii_begin,ii_end
	906	ij=(j-1)*iim+i
	907	phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
	908	ENDDO
	909	ENDDO
	910	ENDIF
	911	!!-> implicit flush on phi(:,1)
	912
	913	!$OMP BARRIER
	914
	915	! for other layers
	916	DO l = ll_beginp1, ll_end
	917	DO j=jj_begin,jj_end
	918	DO i=ii_begin,ii_end
	919	ij=(j-1)*iim+i
	920	phi(ij,l) = 0.5 * ( theta(ij,l) + theta(ij,l-1) ) &
	921	* ( pk(ij,l-1) - pk(ij,l) )
	922	ENDDO
	923	ENDDO
	924	ENDDO
	925
	926	!$OMP BARRIER
	927
	928
	929	IF (is_omp_first_level) THEN
	930	DO l = 2, llm
	931	DO j=jj_begin,jj_end
	932	! ---> Bug compilo intel ici en openmp
	933	! ---> Couper la boucle
	934	IF (j==jj_end+1) PRINT*,"this message must not be printed"
	935	DO i=ii_begin,ii_end
	936	ij=(j-1)*iim+i
	937	phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
	938	ENDDO
	939	ENDDO
	940	ENDDO
	941	! --> IMPLICIT FLUSH on phi --> non
	942	ENDIF
	943
	944	! compute wind centered lon lat compound
	945	DO l=ll_begin,ll_end
	946	DO j=jj_begin,jj_end
	947	DO i=ii_begin,ii_end
	948	ij=(j-1)*iim+i
	949	uc(:)=1/Ai(ij)* &
	950	( 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,:)) &
	951	+ 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,:)) &
	952	+ 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,:)) &
	953	+ 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,:)) &
	954	+ 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,:))&
	955	+ 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,:)))
	956	ulon(ij,l)=sum(uc(:)*elon_i(ij,:))
	957	ulat(ij,l)=sum(uc(:)*elat_i(ij,:))
	958	ENDDO
	959	ENDDO
	960	ENDDO
	961
	962
	963	! compute centered vorticity
	964
	965	DO l=ll_begin,ll_end
	966	DO j=jj_begin,jj_end
	967	DO i=ii_begin,ii_end
	968	ij=(j-1)*iim+i
	969	vortc(ij,l) = Riv(ij,vup) * vort(ij+z_up,l) + &
	970	Riv(ij,vlup) * vort(ij+z_lup,l) + &
	971	Riv(ij,vldown)* vort(ij+z_ldown,l) + &
	972	Riv(ij,vdown) * vort(ij+z_down,l) + &
	973	Riv(ij,vrdown)* vort(ij+z_rdown,l) + &
	974	Riv(ij,vrup) * vort(ij+z_rup,l)
	975	ENDDO
	976	ENDDO
	977	ENDDO
	978
	979
	980
	981	!$OMP BARRIER
	982	END SUBROUTINE grid_icosa_to_physics
	983
	984
	985	SUBROUTINE grid_physics_to_icosa
	986	USE theta2theta_rhodz_mod
	987	USE omp_para
	988	IMPLICIT NONE
	989	INTEGER :: i,j,ij,l,iq
	990
	991	DO l=ll_begin,ll_end
	992	DO j=jj_begin,jj_end
	993	DO i=ii_begin,ii_end
	994	ij=(j-1)*iim+i
	995	duc(ij,:,l)=dulon(ij,l)elon_i(ij,:)+dulat(ij,l)elat_i(ij,:)
	996	ENDDO
	997	ENDDO
	998	ENDDO
	999
	1000	DO l=ll_begin,ll_end
	1001	DO j=jj_begin,jj_end
	1002	DO i=ii_begin,ii_end
	1003	ij=(j-1)*iim+i
	1004	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,:) )
	1005	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,:) )
	1006	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,:) )
	1007	ENDDO
	1008	ENDDO
	1009	ENDDO
	1010
	1011	DO l=ll_begin,ll_end
	1012	DO j=jj_begin,jj_end
	1013	DO i=ii_begin,ii_end
	1014	ij=(j-1)*iim+i
	1015	Temp(ij,l)=Temp(ij,l)+ dtphy * dTemp(ij,l)
	1016	ENDDO
	1017	ENDDO
	1018	ENDDO
	1019
	1020	DO iq=1,nqtot
	1021	DO l=ll_begin,ll_end
	1022	DO j=jj_begin,jj_end
	1023	DO i=ii_begin,ii_end
	1024	ij=(j-1)*iim+i
	1025	q(ij,l,iq)=q(ij,l,iq)+ dtphy * dq(ij,l,iq)
	1026	ENDDO
	1027	ENDDO
	1028	ENDDO
	1029	ENDDO
	1030
	1031	!$OMP BARRIER
	1032
	1033	IF (is_omp_first_level) THEN
	1034	DO j=jj_begin,jj_end
	1035	DO i=ii_begin,ii_end
	1036	ij=(j-1)*iim+i
	1037	ps(ij)=ps(ij)+ dtphy * dps(ij)
	1038	ENDDO
	1039	ENDDO
	1040	ENDIF
	1041
	1042	! CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
	1043
	1044	! compute pression
	1045	!$OMP BARRIER
	1046	DO l = ll_begin,ll_endp1
	1047	DO j=jj_begin,jj_end
	1048	DO i=ii_begin,ii_end
	1049	ij=(j-1)*iim+i
	1050	p(ij,l) = ap(l) + bp(l) * ps(ij)
	1051	ENDDO
	1052	ENDDO
	1053	ENDDO
	1054
	1055	!$OMP BARRIER
	1056
	1057	! compute exner
	1058
	1059	IF (is_omp_first_level) THEN
	1060	DO j=jj_begin,jj_end
	1061	DO i=ii_begin,ii_end
	1062	ij=(j-1)*iim+i
	1063	pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
	1064	ENDDO
	1065	ENDDO
	1066	ENDIF
	1067
	1068	! 3D : pk
	1069	DO l = ll_begin,ll_end
	1070	DO j=jj_begin,jj_end
	1071	DO i=ii_begin,ii_end
	1072	ij=(j-1)*iim+i
	1073	pk(ij,l) = cpp * ((.5/preff)(p(ij,l)+p(ij,l+1))) * kappa
	1074	ENDDO
	1075	ENDDO
	1076	ENDDO
	1077
	1078	!$OMP BARRIER
	1079
	1080	! compute theta, temperature and pression at layer
	1081	DO l = ll_begin, ll_end
	1082	DO j=jj_begin,jj_end
	1083	DO i=ii_begin,ii_end
	1084	ij=(j-1)*iim+i
	1085	theta_rhodz(ij,l,1) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
	1086	ENDDO
	1087	ENDDO
	1088	ENDDO
	1089
	1090	END SUBROUTINE grid_physics_to_icosa
	1091
	1092
	1093
	1094	END SUBROUTINE physics
	1095
	1096
	1097
	1098
	1099
	1100	END MODULE interface_icosa_lmdz_mod

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

Download in other formats:

Original Format