Context Navigation

source: codes/icosagcm/trunk/src/caldyn_gcm.f90 @ 360

Last change on this file since 360 was 360, checked in by dubos, 9 years ago
Towards HEVI time stepping
File size: 34.0 KB

Rev	Line
[12]	1	MODULE caldyn_gcm_mod
[19]	2	USE icosa
[151]	3	USE transfert_mod
[350]	4	IMPLICIT NONE
[132]	5	PRIVATE
	6
	7	INTEGER, PARAMETER :: energy=1, enstrophy=2
[162]	8	TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:)
[186]	9	REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:)
	10	!$OMP THREADPRIVATE(out_u, p, qu)
[17]	11
[151]	12	! temporary shared variable for caldyn
	13	TYPE(t_field),POINTER :: f_pk(:)
	14	TYPE(t_field),POINTER :: f_wwuu(:)
[179]	15	TYPE(t_field),POINTER :: f_planetvel(:)
[151]	16
[159]	17	INTEGER :: caldyn_conserv
[186]	18	!$OMP THREADPRIVATE(caldyn_conserv)
	19
[162]	20	TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
[122]	21
[354]	22	PUBLIC init_caldyn, caldyn_BC, caldyn, req_ps, req_mass
[151]	23
[12]	24	CONTAINS
[15]	25
[98]	26	SUBROUTINE init_caldyn
[50]	27	USE icosa
[354]	28	USE observable_mod
	29	! USE exner_mod
[131]	30	USE mpipara
[295]	31	USE omp_para
[50]	32	IMPLICIT NONE
[122]	33	CHARACTER(len=255) :: def
[179]	34	INTEGER :: ind
	35	REAL(rstd),POINTER :: planetvel(:)
[122]	36
[195]	37	def='energy'
[125]	38	CALL getin('caldyn_conserv',def)
	39	SELECT CASE(TRIM(def))
	40	CASE('energy')
[132]	41	caldyn_conserv=energy
[126]	42	CASE('enstrophy')
[132]	43	caldyn_conserv=enstrophy
[125]	44	CASE DEFAULT
[159]	45	IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
	46	TRIM(def),'> options are <energy>, <enstrophy>'
[125]	47	STOP
	48	END SELECT
[295]	49	IF (is_master) PRINT *, 'caldyn_conserv=',def
[125]	50
[17]	51	CALL allocate_caldyn
[179]	52
	53	DO ind=1,ndomain
[202]	54	IF (.NOT. assigned_domain(ind)) CYCLE
[179]	55	CALL swap_dimensions(ind)
	56	CALL swap_geometry(ind)
	57	planetvel=f_planetvel(ind)
	58	CALL compute_planetvel(planetvel)
	59	END DO
	60
[15]	61	END SUBROUTINE init_caldyn
[179]	62
[12]	63	SUBROUTINE allocate_caldyn
[19]	64	USE icosa
[12]	65	IMPLICIT NONE
	66
[151]	67	CALL allocate_field(f_out_u,field_u,type_real,llm)
	68	CALL allocate_field(f_qu,field_u,type_real,llm)
[162]	69	CALL allocate_field(f_qv,field_z,type_real,llm)
[159]	70	CALL allocate_field(f_pk, field_t,type_real,llm, name='pk')
	71	CALL allocate_field(f_wwuu, field_u,type_real,llm+1,name='wwuu')
[179]	72	CALL allocate_field(f_planetvel, field_u,type_real, name='planetvel') ! planetary velocity at r=a
[151]	73
[12]	74	END SUBROUTINE allocate_caldyn
[157]	75
[350]	76	SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux)
[157]	77	USE icosa
	78	USE mpipara
	79	USE omp_para
	80	TYPE(t_field),POINTER :: f_phis(:)
[350]	81	TYPE(t_field),POINTER :: f_geopot(:)
[157]	82	TYPE(t_field),POINTER :: f_wflux(:)
	83	REAL(rstd),POINTER :: phis(:)
	84	REAL(rstd),POINTER :: wflux(:,:)
	85	REAL(rstd),POINTER :: geopot(:,:)
	86	REAL(rstd),POINTER :: wwuu(:,:)
	87
	88	INTEGER :: ind,i,j,ij,l
	89
[295]	90	IF (is_omp_first_level) THEN
[157]	91	DO ind=1,ndomain
[186]	92	IF (.NOT. assigned_domain(ind)) CYCLE
[157]	93	CALL swap_dimensions(ind)
	94	CALL swap_geometry(ind)
	95	geopot=f_geopot(ind)
	96	phis=f_phis(ind)
	97	wflux=f_wflux(ind)
	98	wwuu=f_wwuu(ind)
	99
[174]	100	DO ij=ij_begin_ext,ij_end_ext
	101	! lower BCs : geopot=phis, wflux=0, wwuu=0
	102	geopot(ij,1) = phis(ij)
	103	wflux(ij,1) = 0.
	104	wwuu(ij+u_right,1)=0
	105	wwuu(ij+u_lup,1)=0
	106	wwuu(ij+u_ldown,1)=0
	107	! top BCs : wflux=0, wwuu=0
	108	wflux(ij,llm+1) = 0.
	109	wwuu(ij+u_right,llm+1)=0
	110	wwuu(ij+u_lup,llm+1)=0
	111	wwuu(ij+u_ldown,llm+1)=0
[157]	112	ENDDO
	113	END DO
	114	ENDIF
	115
[295]	116	!$OMP BARRIER
[157]	117	END SUBROUTINE caldyn_BC
[56]	118
[159]	119	SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
[350]	120	f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
[126]	121	USE icosa
[354]	122	USE observable_mod
[167]	123	USE disvert_mod, ONLY : caldyn_eta, eta_mass
[126]	124	USE vorticity_mod
	125	USE kinetic_mod
	126	USE theta2theta_rhodz_mod
[191]	127	USE wind_mod
[131]	128	USE mpipara
[145]	129	USE trace
[151]	130	USE omp_para
[171]	131	USE output_field_mod
[295]	132	USE checksum_mod
[126]	133	IMPLICIT NONE
[129]	134	LOGICAL,INTENT(IN) :: write_out
[126]	135	TYPE(t_field),POINTER :: f_phis(:)
[12]	136	TYPE(t_field),POINTER :: f_ps(:)
[159]	137	TYPE(t_field),POINTER :: f_mass(:)
[126]	138	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	139	TYPE(t_field),POINTER :: f_u(:)
	140	TYPE(t_field),POINTER :: f_q(:)
[350]	141	TYPE(t_field),POINTER :: f_geopot(:)
[134]	142	TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
[360]	143	TYPE(t_field) :: f_dps(:)
	144	TYPE(t_field) :: f_dmass(:)
	145	TYPE(t_field) :: f_dtheta_rhodz(:)
	146	TYPE(t_field) :: f_du(:)
[12]	147
[159]	148	REAL(rstd),POINTER :: ps(:), dps(:)
	149	REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
[134]	150	REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
[159]	151	REAL(rstd),POINTER :: qu(:,:)
[162]	152	REAL(rstd),POINTER :: qv(:,:)
[151]	153
	154	! temporary shared variable
	155	REAL(rstd),POINTER :: theta(:,:)
[157]	156	REAL(rstd),POINTER :: pk(:,:)
[156]	157	REAL(rstd),POINTER :: geopot(:,:)
[162]	158	REAL(rstd),POINTER :: convm(:,:)
[151]	159	REAL(rstd),POINTER :: wwuu(:,:)
[157]	160
	161	INTEGER :: ind
[151]	162	LOGICAL,SAVE :: first=.TRUE.
	163	!$OMP THREADPRIVATE(first)
[12]	164
[162]	165	! MPI messages need to be sent at first call to caldyn
	166	! This is needed only once : the next ones will be sent by timeloop
	167	IF (first) THEN
[151]	168	first=.FALSE.
[162]	169	IF(caldyn_eta==eta_mass) THEN
	170	CALL init_message(f_ps,req_i1,req_ps)
[190]	171	ELSE
[162]	172	CALL init_message(f_mass,req_i1,req_mass)
	173	END IF
[151]	174	CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
	175	CALL init_message(f_u,req_e1_vect,req_u)
	176	CALL init_message(f_qu,req_e1_scal,req_qu)
[186]	177	! IF(caldyn_eta==eta_mass) THEN
	178	! CALL send_message(f_ps,req_ps)
	179	! CALL wait_message(req_ps)
	180	! ELSE
	181	! CALL send_message(f_mass,req_mass)
	182	! CALL wait_message(req_mass)
	183	! END IF
	184	ENDIF
	185
	186	CALL trace_start("caldyn")
	187
[162]	188	IF(caldyn_eta==eta_mass) THEN
	189	CALL send_message(f_ps,req_ps)
	190	ELSE
	191	CALL send_message(f_mass,req_mass)
	192	END IF
[126]	193
[327]	194	CALL send_message(f_theta_rhodz,req_theta_rhodz)
[151]	195	CALL send_message(f_u,req_u)
	196
[126]	197	SELECT CASE(caldyn_conserv)
[132]	198	CASE(energy) ! energy-conserving
[128]	199	DO ind=1,ndomain
[186]	200	IF (.NOT. assigned_domain(ind)) CYCLE
[128]	201	CALL swap_dimensions(ind)
	202	CALL swap_geometry(ind)
	203	ps=f_ps(ind)
[157]	204	u=f_u(ind)
	205	theta_rhodz = f_theta_rhodz(ind)
[159]	206	mass=f_mass(ind)
[157]	207	theta = f_theta(ind)
[128]	208	qu=f_qu(ind)
[162]	209	qv=f_qv(ind)
	210	CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
[128]	211	ENDDO
[347]	212	! CALL checksum(f_mass)
	213	! CALL checksum(f_theta)
[128]	214
[151]	215	CALL send_message(f_qu,req_qu)
[327]	216	! CALL wait_message(req_qu)
[128]	217
	218	DO ind=1,ndomain
[186]	219	IF (.NOT. assigned_domain(ind)) CYCLE
[128]	220	CALL swap_dimensions(ind)
	221	CALL swap_geometry(ind)
	222	ps=f_ps(ind)
[157]	223	u=f_u(ind)
[128]	224	theta_rhodz=f_theta_rhodz(ind)
[159]	225	mass=f_mass(ind)
[157]	226	theta = f_theta(ind)
[128]	227	qu=f_qu(ind)
[151]	228	pk = f_pk(ind)
[156]	229	geopot = f_geopot(ind)
[183]	230	CALL compute_geopot(ps,mass,theta, pk,geopot)
[157]	231	hflux=f_hflux(ind)
[162]	232	convm = f_dmass(ind)
[157]	233	dtheta_rhodz=f_dtheta_rhodz(ind)
	234	du=f_du(ind)
[162]	235	CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
	236	IF(caldyn_eta==eta_mass) THEN
	237	wflux=f_wflux(ind)
	238	wwuu=f_wwuu(ind)
	239	dps=f_dps(ind)
	240	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
	241	END IF
[128]	242	ENDDO
[347]	243
	244	! CALL checksum(f_geopot)
	245	! CALL checksum(f_dmass)
	246	! CALL checksum(f_pk)
	247	! CALL checksum(f_pk)
[128]	248
[132]	249	CASE(enstrophy) ! enstrophy-conserving
[126]	250	DO ind=1,ndomain
[186]	251	IF (.NOT. assigned_domain(ind)) CYCLE
[126]	252	CALL swap_dimensions(ind)
	253	CALL swap_geometry(ind)
[157]	254	ps=f_ps(ind)
	255	u=f_u(ind)
	256	theta_rhodz=f_theta_rhodz(ind)
[159]	257	mass=f_mass(ind)
[157]	258	theta = f_theta(ind)
	259	qu=f_qu(ind)
[182]	260	qv=f_qv(ind)
[162]	261	CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
[157]	262	pk = f_pk(ind)
[156]	263	geopot = f_geopot(ind)
[183]	264	CALL compute_geopot(ps,mass,theta, pk,geopot)
[134]	265	hflux=f_hflux(ind)
[162]	266	convm = f_dmass(ind)
[126]	267	dtheta_rhodz=f_dtheta_rhodz(ind)
	268	du=f_du(ind)
[162]	269	CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
	270	IF(caldyn_eta==eta_mass) THEN
	271	wflux=f_wflux(ind)
	272	wwuu=f_wwuu(ind)
	273	dps=f_dps(ind)
	274	CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
	275	END IF
[126]	276	ENDDO
	277
	278	CASE DEFAULT
	279	STOP
	280	END SELECT
[12]	281
[295]	282	!$OMP BARRIER
[126]	283	! CALL check_mass_conservation(f_ps,f_dps)
[145]	284	CALL trace_end("caldyn")
[186]	285	!!$OMP BARRIER
[126]	286
	287	END SUBROUTINE caldyn
[179]	288
	289	SUBROUTINE compute_planetvel(planetvel)
	290	USE wind_mod
	291	REAL(rstd),INTENT(OUT) :: planetvel(iim3jjm)
	292	REAL(rstd) :: ulon(iim3jjm)
	293	REAL(rstd) :: ulat(iim3jjm)
	294	REAL(rstd) :: lon,lat
	295	INTEGER :: ij
	296	DO ij=ij_begin_ext,ij_end_ext
[350]	297	ulon(ij+u_right)=radiusomegacos(lat_e(ij+u_right))
[179]	298	ulat(ij+u_right)=0
	299
[350]	300	ulon(ij+u_lup)=radiusomegacos(lat_e(ij+u_lup))
[179]	301	ulat(ij+u_lup)=0
	302
[350]	303	ulon(ij+u_ldown)=radiusomegacos(lat_e(ij+u_ldown))
[179]	304	ulat(ij+u_ldown)=0
	305	END DO
	306	CALL compute_wind2D_perp_from_lonlat_compound(ulon, ulat, planetvel)
	307	END SUBROUTINE compute_planetvel
[128]	308
[162]	309	SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
[19]	310	USE icosa
[167]	311	USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
[50]	312	USE exner_mod
[145]	313	USE trace
[151]	314	USE omp_para
[12]	315	IMPLICIT NONE
[128]	316	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
	317	REAL(rstd),INTENT(IN) :: ps(iim*jjm)
[157]	318	REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
[162]	319	REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
[157]	320	REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
[128]	321	REAL(rstd),INTENT(OUT) :: qu(iim3jjm,llm)
[162]	322	REAL(rstd),INTENT(OUT) :: qv(iim2jjm,llm)
[128]	323
	324	INTEGER :: i,j,ij,l
[162]	325	REAL(rstd) :: etav,hv, m
	326	! REAL(rstd) :: qv(2iimjjm,llm) ! potential velocity
[128]	327
[151]	328	CALL trace_start("compute_pvort")
[145]	329
[162]	330	IF(caldyn_eta==eta_mass) THEN
[327]	331	CALL wait_message(req_ps)
[162]	332	ELSE
[327]	333	CALL wait_message(req_mass)
[162]	334	END IF
[327]	335	CALL wait_message(req_theta_rhodz)
[151]	336
[162]	337	IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
	338	DO l = ll_begin,ll_end
[327]	339	CALL test_message(req_u)
[174]	340	!DIR$ SIMD
	341	DO ij=ij_begin_ext,ij_end_ext
	342	m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
	343	rhodz(ij,l) = m
	344	theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
[128]	345	ENDDO
	346	ENDDO
[162]	347	ELSE ! Compute only theta
	348	DO l = ll_begin,ll_end
[327]	349	CALL test_message(req_u)
[174]	350	!DIR$ SIMD
	351	DO ij=ij_begin_ext,ij_end_ext
	352	theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
	353	ENDDO
[162]	354	ENDDO
	355	END IF
[151]	356
[327]	357	CALL wait_message(req_u)
[128]	358
[123]	359	!!! Compute shallow-water potential vorticity
[151]	360	DO l = ll_begin,ll_end
[174]	361	!DIR$ SIMD
[179]	362	DO ij=ij_begin_ext,ij_end_ext
[151]	363	etav= 1./Av(ij+z_up)( ne_rup u(ij+u_rup,l) * de(ij+u_rup) &
	364	+ ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left) &
	365	- ne_lup * u(ij+u_lup,l) * de(ij+u_lup) )
[123]	366
	367	hv = Riv2(ij,vup) * rhodz(ij,l) &
	368	+ Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) &
	369	+ Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
	370
	371	qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
	372
[151]	373	etav = 1./Av(ij+z_down)( ne_ldown u(ij+u_ldown,l) * de(ij+u_ldown) &
	374	+ ne_right * u(ij+t_ldown+u_right,l) * de(ij+t_ldown+u_right) &
	375	- ne_rdown * u(ij+u_rdown,l) * de(ij+u_rdown) )
[123]	376
	377	hv = Riv2(ij,vdown) * rhodz(ij,l) &
	378	+ Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) &
	379	+ Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
	380
	381	qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
	382
[12]	383	ENDDO
	384
[174]	385	!DIR$ SIMD
	386	DO ij=ij_begin,ij_end
	387	qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))
	388	qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))
	389	qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))
[126]	390	END DO
	391
	392	ENDDO
	393
[151]	394	CALL trace_end("compute_pvort")
[145]	395
[126]	396	END SUBROUTINE compute_pvort
[125]	397
[183]	398	SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot)
[126]	399	USE icosa
	400	USE disvert_mod
	401	USE exner_mod
[145]	402	USE trace
[151]	403	USE omp_para
[126]	404	IMPLICIT NONE
[159]	405	REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
	406	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
	407	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) ! potential temperature
	408	REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
	409	REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
[157]	410
	411	INTEGER :: i,j,ij,l
	412	REAL(rstd) :: p_ik, exner_ik
[347]	413	INTEGER :: ij_omp_begin_ext, ij_omp_end_ext
[157]	414
[327]	415
[157]	416	CALL trace_start("compute_geopot")
[327]	417
[347]	418	CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
	419	ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
	420	ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
[157]	421
[159]	422	IF(caldyn_eta==eta_mass) THEN
[181]	423
[157]	424	!!! Compute exner function and geopotential
[295]	425	DO l = 1,llm
[181]	426	!DIR$ SIMD
[327]	427	DO ij=ij_omp_begin_ext,ij_omp_end_ext
[295]	428	p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
	429	! p_ik = ptop + g(mass_ak(l)+ mass_bk(l)ps(i,j))
	430	exner_ik = cpp * (p_ik/preff) ** kappa
	431	pk(ij,l) = exner_ik
[181]	432	! specific volume v = kappathetapi/p = dphi/g/rhodz
[295]	433	geopot(ij,l+1) = geopot(ij,l) + (gkappa)rhodz(ij,l)theta(ij,l)exner_ik/p_ik
[159]	434	ENDDO
[295]	435	ENDDO
[327]	436	! ENDIF
[167]	437	ELSE
[159]	438	! We are using a Lagrangian vertical coordinate
	439	! Pressure must be computed first top-down (temporarily stored in pk)
	440	! Then Exner pressure and geopotential are computed bottom-up
[167]	441	! Notice that the computation below should work also when caldyn_eta=eta_mass
[159]	442
[183]	443	IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz
[181]	444	! specific volume 1 = dphi/g/rhodz
[327]	445	! IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
	446	DO l = 1,llm
	447	!DIR$ SIMD
	448	DO ij=ij_omp_begin_ext,ij_omp_end_ext
	449	geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
[295]	450	ENDDO
[327]	451	ENDDO
[183]	452	ELSE ! non-Boussinesq, compute geopotential and Exner pressure
[181]	453	! uppermost layer
[186]	454
[327]	455	!DIR$ SIMD
	456	DO ij=ij_omp_begin_ext,ij_omp_end_ext
	457	pk(ij,llm) = ptop + (.5g)rhodz(ij,llm)
	458	END DO
	459	! other layers
	460	DO l = llm-1, 1, -1
	461	!DIR$ SIMD
	462	DO ij=ij_omp_begin_ext,ij_omp_end_ext
	463	pk(ij,l) = pk(ij,l+1) + (.5g)(rhodz(ij,l)+rhodz(ij,l+1))
	464	END DO
	465	END DO
	466	! surface pressure (for diagnostics)
	467	DO ij=ij_omp_begin_ext,ij_omp_end_ext
	468	ps(ij) = pk(ij,1) + (.5g)rhodz(ij,1)
	469	END DO
	470
	471	! specific volume v = kappathetapi/p = dphi/g/rhodz
	472	DO l = 1,llm
[295]	473	!DIR$ SIMD
[327]	474	DO ij=ij_omp_begin_ext,ij_omp_end_ext
	475	p_ik = pk(ij,l)
	476	exner_ik = cpp * (p_ik/preff) ** kappa
	477	geopot(ij,l+1) = geopot(ij,l) + (gkappa)rhodz(ij,l)theta(ij,l)exner_ik/p_ik
	478	pk(ij,l) = exner_ik
[295]	479	ENDDO
[327]	480	ENDDO
[167]	481	END IF
[157]	482
[159]	483	END IF
	484
[295]	485	!ym flush geopot
	486	!$OMP BARRIER
	487
[157]	488	CALL trace_end("compute_geopot")
	489
	490	END SUBROUTINE compute_geopot
	491
[162]	492	SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du)
[157]	493	USE icosa
	494	USE disvert_mod
	495	USE exner_mod
	496	USE trace
	497	USE omp_para
	498	IMPLICIT NONE
[183]	499	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm) ! prognostic "velocity"
[157]	500	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[126]	501	REAL(rstd),INTENT(IN) :: qu(iim3jjm,llm)
[157]	502	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm) ! potential temperature
[167]	503	REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
[157]	504	REAL(rstd),INTENT(IN) :: geopot(iim*jjm,llm+1) ! geopotential
[126]	505
[157]	506	REAL(rstd),INTENT(OUT) :: hflux(iim3jjm,llm) ! hflux in kg/s
[162]	507	REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm) ! mass flux convergence
[126]	508	REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
[157]	509	REAL(rstd),INTENT(OUT) :: du(iim3jjm,llm)
[151]	510
[183]	511	REAL(rstd) :: cor_NT(iim*jjm,llm) ! NT coriolis force u.(du/dPhi)
[179]	512	REAL(rstd) :: urel(3iimjjm,llm) ! relative velocity
[151]	513	REAL(rstd) :: Ftheta(3iimjjm,llm) ! theta flux
[157]	514	REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function
[126]	515
[139]	516	INTEGER :: i,j,ij,l
[157]	517	REAL(rstd) :: ww,uu
[139]	518
[157]	519	CALL trace_start("compute_caldyn_horiz")
[126]	520
[186]	521	! CALL wait_message(req_theta_rhodz)
[151]	522
	523	DO l = ll_begin, ll_end
[123]	524	!!! Compute mass and theta fluxes
[327]	525	IF (caldyn_conserv==energy) CALL test_message(req_qu)
[174]	526	!DIR$ SIMD
	527	DO ij=ij_begin_ext,ij_end_ext
[134]	528	hflux(ij+u_right,l)=0.5(rhodz(ij,l)+rhodz(ij+t_right,l))u(ij+u_right,l)*le(ij+u_right)
	529	hflux(ij+u_lup,l)=0.5(rhodz(ij,l)+rhodz(ij+t_lup,l))u(ij+u_lup,l)*le(ij+u_lup)
	530	hflux(ij+u_ldown,l)=0.5(rhodz(ij,l)+rhodz(ij+t_ldown,l))u(ij+u_ldown,l)*le(ij+u_ldown)
[12]	531
[134]	532	Ftheta(ij+u_right,l)=0.5(theta(ij,l)+theta(ij+t_right,l))hflux(ij+u_right,l)
	533	Ftheta(ij+u_lup,l)=0.5(theta(ij,l)+theta(ij+t_lup,l))hflux(ij+u_lup,l)
	534	Ftheta(ij+u_ldown,l)=0.5(theta(ij,l)+theta(ij+t_ldown,l))hflux(ij+u_ldown,l)
[12]	535	ENDDO
[123]	536
	537	!!! compute horizontal divergence of fluxes
[174]	538	!DIR$ SIMD
	539	DO ij=ij_begin,ij_end
[162]	540	! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
	541	convm(ij,l)= -1./Ai(ij)(ne_righthflux(ij+u_right,l) + &
[151]	542	ne_rup*hflux(ij+u_rup,l) + &
	543	ne_lup*hflux(ij+u_lup,l) + &
	544	ne_left*hflux(ij+u_left,l) + &
	545	ne_ldown*hflux(ij+u_ldown,l) + &
	546	ne_rdown*hflux(ij+u_rdown,l))
[123]	547
	548	! signe ? attention d (rho theta dz)
[22]	549	! dtheta_rhodz = -div(flux.theta)
[151]	550	dtheta_rhodz(ij,l)=-1./Ai(ij)(ne_rightFtheta(ij+u_right,l) + &
	551	ne_rup*Ftheta(ij+u_rup,l) + &
	552	ne_lup*Ftheta(ij+u_lup,l) + &
	553	ne_left*Ftheta(ij+u_left,l) + &
	554	ne_ldown*Ftheta(ij+u_ldown,l) + &
	555	ne_rdown*Ftheta(ij+u_rdown,l))
[12]	556	ENDDO
[151]	557
[157]	558	END DO
[151]	559
[56]	560	!!! Compute potential vorticity (Coriolis) contribution to du
[157]	561
[128]	562	SELECT CASE(caldyn_conserv)
[132]	563	CASE(energy) ! energy-conserving TRiSK
[12]	564
[327]	565	CALL wait_message(req_qu)
[151]	566
	567	DO l=ll_begin,ll_end
[174]	568	!DIR$ SIMD
	569	DO ij=ij_begin,ij_end
	570
[134]	571	uu = wee(ij+u_right,1,1)hflux(ij+u_rup,l)(qu(ij+u_right,l)+qu(ij+u_rup,l))+ &
	572	wee(ij+u_right,2,1)hflux(ij+u_lup,l)(qu(ij+u_right,l)+qu(ij+u_lup,l))+ &
	573	wee(ij+u_right,3,1)hflux(ij+u_left,l)(qu(ij+u_right,l)+qu(ij+u_left,l))+ &
	574	wee(ij+u_right,4,1)hflux(ij+u_ldown,l)(qu(ij+u_right,l)+qu(ij+u_ldown,l))+ &
	575	wee(ij+u_right,5,1)hflux(ij+u_rdown,l)(qu(ij+u_right,l)+qu(ij+u_rdown,l))+ &
	576	wee(ij+u_right,1,2)hflux(ij+t_right+u_ldown,l)(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+ &
	577	wee(ij+u_right,2,2)hflux(ij+t_right+u_rdown,l)(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+ &
	578	wee(ij+u_right,3,2)hflux(ij+t_right+u_right,l)(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+ &
	579	wee(ij+u_right,4,2)hflux(ij+t_right+u_rup,l)(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+ &
	580	wee(ij+u_right,5,2)hflux(ij+t_right+u_lup,l)(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
[128]	581	du(ij+u_right,l) = .5*uu/de(ij+u_right)
	582
[134]	583	uu = wee(ij+u_lup,1,1)hflux(ij+u_left,l)(qu(ij+u_lup,l)+qu(ij+u_left,l)) + &
	584	wee(ij+u_lup,2,1)hflux(ij+u_ldown,l)(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) + &
	585	wee(ij+u_lup,3,1)hflux(ij+u_rdown,l)(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) + &
	586	wee(ij+u_lup,4,1)hflux(ij+u_right,l)(qu(ij+u_lup,l)+qu(ij+u_right,l)) + &
	587	wee(ij+u_lup,5,1)hflux(ij+u_rup,l)(qu(ij+u_lup,l)+qu(ij+u_rup,l)) + &
	588	wee(ij+u_lup,1,2)hflux(ij+t_lup+u_right,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) + &
	589	wee(ij+u_lup,2,2)hflux(ij+t_lup+u_rup,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) + &
	590	wee(ij+u_lup,3,2)hflux(ij+t_lup+u_lup,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + &
	591	wee(ij+u_lup,4,2)hflux(ij+t_lup+u_left,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) + &
	592	wee(ij+u_lup,5,2)hflux(ij+t_lup+u_ldown,l)(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
[128]	593	du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
[12]	594
[128]	595
[134]	596	uu = wee(ij+u_ldown,1,1)hflux(ij+u_rdown,l)(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) + &
	597	wee(ij+u_ldown,2,1)hflux(ij+u_right,l)(qu(ij+u_ldown,l)+qu(ij+u_right,l)) + &
	598	wee(ij+u_ldown,3,1)hflux(ij+u_rup,l)(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) + &
	599	wee(ij+u_ldown,4,1)hflux(ij+u_lup,l)(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) + &
	600	wee(ij+u_ldown,5,1)hflux(ij+u_left,l)(qu(ij+u_ldown,l)+qu(ij+u_left,l)) + &
	601	wee(ij+u_ldown,1,2)hflux(ij+t_ldown+u_lup,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) + &
	602	wee(ij+u_ldown,2,2)hflux(ij+t_ldown+u_left,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) + &
	603	wee(ij+u_ldown,3,2)hflux(ij+t_ldown+u_ldown,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) + &
	604	wee(ij+u_ldown,4,2)hflux(ij+t_ldown+u_rdown,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) + &
	605	wee(ij+u_ldown,5,2)hflux(ij+t_ldown+u_right,l)(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
[128]	606	du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
	607
	608	ENDDO
	609	ENDDO
[146]	610
[132]	611	CASE(enstrophy) ! enstrophy-conserving TRiSK
[128]	612
[151]	613	DO l=ll_begin,ll_end
[174]	614	!DIR$ SIMD
	615	DO ij=ij_begin,ij_end
[12]	616
[134]	617	uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+ &
	618	wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+ &
	619	wee(ij+u_right,3,1)*hflux(ij+u_left,l)+ &
	620	wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+ &
	621	wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+ &
	622	wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+ &
	623	wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+ &
	624	wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+ &
	625	wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+ &
	626	wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
[128]	627	du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
	628
	629
[134]	630	uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+ &
	631	wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+ &
	632	wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+ &
	633	wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+ &
	634	wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+ &
	635	wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+ &
	636	wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+ &
	637	wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+ &
	638	wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+ &
	639	wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
[128]	640	du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
	641
[134]	642	uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+ &
	643	wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+ &
	644	wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+ &
	645	wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+ &
	646	wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+ &
	647	wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+ &
	648	wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+ &
	649	wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+ &
	650	wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+ &
	651	wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
[128]	652	du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
[12]	653
[128]	654	ENDDO
	655	ENDDO
	656
	657	CASE DEFAULT
	658	STOP
	659	END SELECT
[12]	660
	661	!!! Compute bernouilli term = Kinetic Energy + geopotential
[167]	662	IF(boussinesq) THEN
[183]	663	! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)
	664	! uppermost layer
	665	!DIR$ SIMD
	666	DO ij=ij_begin_ext,ij_end_ext
	667	pk(ij,llm) = ptop + (.5g)theta(ij,llm)*rhodz(ij,llm)
	668	END DO
	669	! other layers
	670	DO l = llm-1, 1, -1
[186]	671	! !$OMP DO SCHEDULE(STATIC)
[183]	672	!DIR$ SIMD
	673	DO ij=ij_begin_ext,ij_end_ext
	674	pk(ij,l) = pk(ij,l+1) + (.5g)(theta(ij,l)rhodz(ij,l)+theta(ij,l+1)rhodz(ij,l+1))
	675	END DO
	676	END DO
	677	! surface pressure (for diagnostics) FIXME
	678	! DO ij=ij_begin_ext,ij_end_ext
	679	! ps(ij) = pk(ij,1) + (.5g)theta(ij,1)*rhodz(ij,1)
	680	! END DO
	681	! now pk contains the Lagrange multiplier (pressure)
[12]	682
[167]	683	DO l=ll_begin,ll_end
[174]	684	!DIR$ SIMD
	685	DO ij=ij_begin,ij_end
[186]	686
	687	berni(ij,l) = pk(ij,l) + &
	688	1/(4Ai(ij))(le(ij+u_right)de(ij+u_right)u(ij+u_right,l)**2 + &
[12]	689	le(ij+u_rup)de(ij+u_rup)u(ij+u_rup,l)**2 + &
	690	le(ij+u_lup)de(ij+u_lup)u(ij+u_lup,l)**2 + &
	691	le(ij+u_left)de(ij+u_left)u(ij+u_left,l)**2 + &
	692	le(ij+u_ldown)de(ij+u_ldown)u(ij+u_ldown,l)**2 + &
	693	le(ij+u_rdown)de(ij+u_rdown)u(ij+u_rdown,l)**2 )
[181]	694	! from now on pk contains the vertically-averaged geopotential
[167]	695	pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
	696	ENDDO
	697	ENDDO
	698
[183]	699	ELSE ! compressible
[167]	700
	701	DO l=ll_begin,ll_end
[174]	702	!DIR$ SIMD
	703	DO ij=ij_begin,ij_end
[167]	704
	705	berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
	706	+ 1/(4Ai(ij))(le(ij+u_right)de(ij+u_right)u(ij+u_right,l)**2 + &
	707	le(ij+u_rup)de(ij+u_rup)u(ij+u_rup,l)**2 + &
	708	le(ij+u_lup)de(ij+u_lup)u(ij+u_lup,l)**2 + &
	709	le(ij+u_left)de(ij+u_left)u(ij+u_left,l)**2 + &
	710	le(ij+u_ldown)de(ij+u_ldown)u(ij+u_ldown,l)**2 + &
	711	le(ij+u_rdown)de(ij+u_rdown)u(ij+u_rdown,l)**2 )
	712	ENDDO
	713	ENDDO
[151]	714
[167]	715	END IF ! Boussinesq/compressible
	716
[157]	717	!!! Add gradients of Bernoulli and Exner functions to du
[167]	718	DO l=ll_begin,ll_end
[174]	719	!DIR$ SIMD
	720	DO ij=ij_begin,ij_end
[12]	721
[167]	722	du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * ( &
	723	0.5*(theta(ij,l)+theta(ij+t_right,l)) &
	724	( ne_rightpk(ij,l)+ne_left*pk(ij+t_right,l)) &
	725	+ ne_rightberni(ij,l)+ne_leftberni(ij+t_right,l) )
[12]	726
[123]	727
[167]	728	du(ij+u_lup,l) = du(ij+u_lup,l) + 1/de(ij+u_lup) * ( &
	729	0.5*(theta(ij,l)+theta(ij+t_lup,l)) &
	730	( ne_luppk(ij,l)+ne_rdown*pk(ij+t_lup,l)) &
	731	+ ne_lupberni(ij,l)+ne_rdownberni(ij+t_lup,l) )
[123]	732
[167]	733	du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * ( &
	734	0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
	735	( ne_ldownpk(ij,l)+ne_rup*pk(ij+t_ldown,l)) &
	736	+ ne_ldownberni(ij,l)+ne_rupberni(ij+t_ldown,l) )
[123]	737
[167]	738	ENDDO
[12]	739	ENDDO
[167]	740
	741	CALL trace_end("compute_caldyn_horiz")
[151]	742
[157]	743	END SUBROUTINE compute_caldyn_horiz
	744
[162]	745	SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
[157]	746	USE icosa
	747	USE disvert_mod
	748	USE exner_mod
	749	USE trace
	750	USE omp_para
	751	IMPLICIT NONE
	752	REAL(rstd),INTENT(IN) :: u(iim3jjm,llm)
	753	REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
	754	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[162]	755	REAL(rstd),INTENT(INOUT) :: convm(iim*jjm,llm) ! mass flux convergence
[157]	756
[186]	757	REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
	758	REAL(rstd),INTENT(INOUT) :: wwuu(iim3jjm,llm+1)
[327]	759	REAL(rstd),INTENT(INOUT) :: du(iim3jjm,llm)
	760	REAL(rstd),INTENT(INOUT) :: dtheta_rhodz(iim*jjm,llm)
[157]	761	REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
	762
	763	! temporary variable
	764	INTEGER :: i,j,ij,l
	765	REAL(rstd) :: p_ik, exner_ik
[347]	766	INTEGER :: ij_omp_begin, ij_omp_end
[157]	767
[327]	768
	769	CALL trace_start("compute_geopot")
	770
[347]	771	CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end)
	772	ij_omp_begin=ij_omp_begin+ij_begin-1
	773	ij_omp_end=ij_omp_end+ij_begin-1
	774
[186]	775	! REAL(rstd) :: wwuu(iim3jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp
	776	! need to be understood
[174]	777
[186]	778	! wwuu=wwuu_out
[157]	779	CALL trace_start("compute_caldyn_vert")
	780
[295]	781	!$OMP BARRIER
[162]	782	!!! cumulate mass flux convergence from top to bottom
[327]	783	! IF (is_omp_level_master) THEN
[295]	784	DO l = llm-1, 1, -1
[186]	785	! IF (caldyn_conserv==energy) CALL test_message(req_qu)
	786
	787	!!$OMP DO SCHEDULE(STATIC)
[174]	788	!DIR$ SIMD
[327]	789	DO ij=ij_omp_begin,ij_omp_end
[295]	790	convm(ij,l) = convm(ij,l) + convm(ij,l+1)
	791	ENDDO
[151]	792	ENDDO
[327]	793	! ENDIF
[157]	794
[295]	795	!$OMP BARRIER
	796	! FLUSH on convm
[157]	797	!!!!!!!!!!!!!!!!!!!!!!!!!
	798
	799	! compute dps
[295]	800	IF (is_omp_first_level) THEN
[174]	801	!DIR$ SIMD
	802	DO ij=ij_begin,ij_end
[157]	803	! dps/dt = -int(div flux)dz
[162]	804	dps(ij) = convm(ij,1) * g
[157]	805	ENDDO
	806	ENDIF
[151]	807
[157]	808	!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
[151]	809	DO l=ll_beginp1,ll_end
[186]	810	! IF (caldyn_conserv==energy) CALL test_message(req_qu)
[174]	811	!DIR$ SIMD
	812	DO ij=ij_begin,ij_end
[157]	813	! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
	814	! => w>0 for upward transport
[162]	815	wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l )
[151]	816	ENDDO
	817	ENDDO
[157]	818
[295]	819	!--> flush wflux
	820	!$OMP BARRIER
	821
[157]	822	DO l=ll_begin,ll_endm1
[174]	823	!DIR$ SIMD
	824	DO ij=ij_begin,ij_end
[157]	825	dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) - 0.5 * ( wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
[151]	826	ENDDO
[157]	827	ENDDO
	828
	829	DO l=ll_beginp1,ll_end
[174]	830	!DIR$ SIMD
	831	DO ij=ij_begin,ij_end
[157]	832	dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l ) + 0.5 * ( wflux(ij,l ) * (theta(ij,l-1) + theta(ij,l) ) )
[151]	833	ENDDO
[157]	834	ENDDO
[295]	835
[151]	836
[157]	837	! Compute vertical transport
[151]	838	DO l=ll_beginp1,ll_end
[174]	839	!DIR$ SIMD
	840	DO ij=ij_begin,ij_end
[151]	841	wwuu(ij+u_right,l) = 0.5( wflux(ij,l) + wflux(ij+t_right,l)) (u(ij+u_right,l) - u(ij+u_right,l-1))
	842	wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
	843	wwuu(ij+u_ldown,l) = 0.5( wflux(ij,l) + wflux(ij+t_ldown,l)) (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
	844	ENDDO
	845	ENDDO
[12]	846
[151]	847	!--> flush wwuu
[295]	848	!$OMP BARRIER
[12]	849
[157]	850	! Add vertical transport to du
[151]	851	DO l=ll_begin,ll_end
[174]	852	!DIR$ SIMD
	853	DO ij=ij_begin,ij_end
[151]	854	du(ij+u_right, l ) = du(ij+u_right,l) - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
	855	du(ij+u_lup, l ) = du(ij+u_lup,l) - (wwuu(ij+u_lup,l+1) + wwuu(ij+u_lup,l)) / (rhodz(ij,l)+rhodz(ij+t_lup,l))
	856	du(ij+u_ldown, l ) = du(ij+u_ldown,l) - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
[12]	857	ENDDO
	858	ENDDO
[186]	859
	860	! DO l=ll_beginp1,ll_end
	861	!!DIR$ SIMD
	862	! DO ij=ij_begin,ij_end
	863	! wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l)
	864	! wwuu_out(ij+u_lup,l) = wwuu(ij+u_lup,l)
	865	! wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l)
	866	! ENDDO
	867	! ENDDO
[151]	868
[157]	869	CALL trace_end("compute_caldyn_vert")
[145]	870
[157]	871	END SUBROUTINE compute_caldyn_vert
[126]	872
	873	!-------------------------------- Diagnostics ----------------------------
	874
	875	SUBROUTINE check_mass_conservation(f_ps,f_dps)
	876	USE icosa
[131]	877	USE mpipara
[126]	878	IMPLICIT NONE
	879	TYPE(t_field),POINTER :: f_ps(:)
	880	TYPE(t_field),POINTER :: f_dps(:)
	881	REAL(rstd),POINTER :: ps(:)
	882	REAL(rstd),POINTER :: dps(:)
	883	REAL(rstd) :: mass_tot,dmass_tot
	884	INTEGER :: ind,i,j,ij
	885
	886	mass_tot=0
	887	dmass_tot=0
	888
	889	CALL transfert_request(f_dps,req_i1)
	890	CALL transfert_request(f_ps,req_i1)
	891
	892	DO ind=1,ndomain
	893	CALL swap_dimensions(ind)
	894	CALL swap_geometry(ind)
	895
	896	ps=f_ps(ind)
	897	dps=f_dps(ind)
	898
	899	DO j=jj_begin,jj_end
	900	DO i=ii_begin,ii_end
	901	ij=(j-1)*iim+i
	902	IF (domain(ind)%own(i,j)) THEN
	903	mass_tot=mass_tot+ps(ij)*Ai(ij)/g
	904	dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
	905	ENDIF
	906	ENDDO
	907	ENDDO
	908
	909	ENDDO
[131]	910	IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot," dmass_tot ",dmass_tot
[126]	911
	912	END SUBROUTINE check_mass_conservation
[12]	913
	914	END MODULE caldyn_gcm_mod

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

Download in other formats: