Context Navigation

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

Last change on this file since 176 was 174, checked in by ymipsl, 11 years ago

Transform 2 loops on i and j in one loop ij for efficiency, vectorization and future GPU programing

YM

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

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

Download in other formats: