Context Navigation

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

Last change on this file since 275 was 266, checked in by ymipsl, 10 years ago

Synchronize trunk and Saturn branch.
Merge modification from Saturn branch to trunk

YM

File size: 39.4 KB

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

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

Download in other formats: