Context Navigation

source: codes/icosagcm/trunk/src/diagnostics/check_conserve.f90 @ 649

Last change on this file since 649 was 649, checked in by aslmd, 7 years ago
added angular momentum diagnostics associated with either positive or negative zonal wind, to diagnose the true impact of possible AAM spurious drift on the resulting circulations
File size: 16.0 KB

Rev	Line
[149]	1	MODULE check_conserve_mod
[151]	2	USE icosa
	3	IMPLICIT NONE
[149]	4
[365]	5	PRIVATE
	6
	7	INTEGER, PARAMETER, PUBLIC :: check_basic=1, check_detailed=2, &
	8	AAM_dissip=1, AAM_dyn=2, AAM_phys=3
	9	INTEGER, SAVE :: check_type
[326]	10
[365]	11	TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
	12	TYPE(t_field),POINTER,SAVE :: f_vort(:)
	13	TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
	14
	15	REAL(rstd),SAVE :: AAM_mass, AAM_mass_old, AAM_vel, AAM_vel_old, &
[649]	16	AAM_velp, AAM_velp_old, AAM_velm, AAM_velm_old, &
[365]	17	AAM_mass_source(3), AAM_vel_source(3) ! read/written only IF is_master
[649]	18	REAL(rstd),SAVE :: AAM_vel_plus_source(3), AAM_vel_minus_source(3)
[365]	19	REAL(rstd),SAVE :: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0
	20	!$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
	21
	22	PUBLIC :: init_check_conserve, check_conserve_detailed, check_conserve
[151]	23	CONTAINS
[149]	24
	25	!---------------------------------------------------------------------
	26
[151]	27	SUBROUTINE init_check_conserve
[326]	28	USE getin_mod
	29	USE omp_para, ONLY : is_master
	30	CHARACTER(LEN=255) :: check_type_str
[149]	31	CALL allocate_field(f_pk,field_t,type_real,llm)
	32	CALL allocate_field(f_p,field_t,type_real,llm+1)
	33	CALL allocate_field(f_pks,field_t,type_real)
	34	CALL allocate_field(f_rhodz,field_t,type_real,llm)
	35	CALL allocate_field(f_vort,field_z,type_real,llm)
[326]	36
	37	check_type_str='basic'
	38	CALL getin("check_conservation",check_type_str)
	39	SELECT CASE(TRIM(check_type_str))
	40	CASE('basic')
	41	check_type=check_basic
	42	CASE('detailed')
	43	check_type=check_detailed
	44	CASE DEFAULT
	45	IF(is_master) PRINT*,'Bad selector <', TRIM(check_type_str), &
	46	'> for variable check_conservation : options are <basic>, <detailed>'
	47	STOP
	48	END SELECT
[151]	49	END SUBROUTINE init_check_conserve
[149]	50
[326]	51	!--------------------------------- Basic check --------------------------------
	52
[151]	53	SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it)
	54	USE pression_mod
	55	USE vorticity_mod
	56	USE caldyn_gcm_mod
	57	USE exner_mod
[365]	58	! USE mpipara, ONLY : is_mpi_master, comm_icosa
	59	USE omp_para, ONLY : is_master
[151]	60	TYPE(t_field),POINTER :: f_ps(:)
	61	TYPE(t_field),POINTER :: f_dps(:)
	62	TYPE(t_field),POINTER :: f_ue(:)
	63	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	64	TYPE(t_field),POINTER :: f_phis(:)
[326]	65	INTEGER, INTENT(IN) :: it
[160]	66
[151]	67	REAL(rstd),POINTER :: p(:,:),rhodz(:,:)
[365]	68	INTEGER :: ind,ierr
[326]	69	REAL(rstd) :: mtot, angtot, rmsdpdt
[649]	70	REAL(rstd) :: etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot, ztot
[151]	71
[295]	72	CALL transfert_request(f_ue,req_e1_vect)
[151]	73	CALL pression(f_ps,f_p)
	74
	75	DO ind=1,ndomain
[186]	76	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	77	CALL swap_dimensions(ind)
	78	CALL swap_geometry(ind)
[160]	79	p=f_p(ind)
	80	rhodz=f_rhodz(ind)
[151]	81	CALL compute_rhodz(p,rhodz)
[160]	82	END DO
[151]	83
	84	CALL vorticity(f_ue,f_vort)
[160]	85	CALL check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
[295]	86	CALL check_PV(ztot)
[151]	87	CALL exner(f_ps,f_p,f_pks,f_pk)
[649]	88	CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot)
[160]	89
[365]	90	IF (is_master) THEN ! is_master = is_omp_master && is_mpi_master
	91	! AAM_mass and AAM_vel must be read/written only IF is_master
	92	AAM_mass = ang_mass
	93	AAM_vel = ang_vel
[649]	94	AAM_velp = ang_velp
	95	AAM_velm = ang_velm
[365]	96	angtot = ang_mass + ang_vel
[295]	97	IF ( it == itau0 ) THEN
[160]	98	ztot0 = ztot
	99	mtot0 = mtot
	100	etot0 = etot
[326]	101	angtot0 = angtot
[160]	102	stot0 = stot
[365]	103	AAM_mass_old=AAM_mass
	104	AAM_vel_old=AAM_vel
[649]	105	AAM_velp_old=AAM_velp
	106	AAM_velm_old=AAM_velm
[365]	107	AAM_mass_source=0.
	108	AAM_vel_source=0.
[649]	109	AAM_vel_plus_source=0.
	110	AAM_vel_minus_source=0.
[160]	111	END IF
[295]	112	rmsvtot=SQRT(rmsvtot/mtot)
[160]	113	ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1.
[295]	114	etot=etot/etot0-1. ; angtot=angtot/angtot0-1. ; stot=stot/stot0-1.
[160]	115	rmsdpdt= daysec1.e-2sqrt(rmsdpdt/ncell_glo)
	116
[151]	117	OPEN(134,file="checkconsicosa.txt",position='append')
[295]	118	WRITE(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
	119	WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
[151]	120	WRITE(134,*)"=================================================="
[295]	121	WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
[160]	122
[186]	123	4000 FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie' &
[365]	124	,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB ' &
[295]	125	,e10.3,e13.6,5e13.3/)
[365]	126	CLOSE(134)
	127
	128	IF(check_type == check_detailed) THEN
	129	!ang_mass = AAM_mass-AAM_mass_old
	130	!WRITE(*,'(A,3E12.4)') 'AAM_mass sanity check', SUM(AAM_mass_source), ang_mass, SUM(AAM_mass_source)-ang_mass
	131	!ang_vel = AAM_vel-AAM_vel_old
	132	!WRITE(*,'(A,3E12.4)') 'AAM_vel sanity check', SUM(AAM_vel_source), ang_vel, SUM(AAM_vel_source)-ang_vel
	133	WRITE(,'(A,6E12.4)') 'AAM_mass : time,old,new,dissip,dyn,phys', dtit, AAM_mass_old, AAM_mass, &
	134	AAM_mass_source(AAM_dissip), AAM_mass_source(AAM_dyn), AAM_mass_source(AAM_phys)
	135	WRITE(,'(A,6E12.4)') 'AAM_vel : time,old,new,dissip,dyn,phys', dtit, AAM_vel_old, AAM_vel, &
	136	AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys)
[649]	137	WRITE(,'(A,6E12.4)') 'AAM_vel+ : time,old,new,dissip,dyn,phys', dtit, AAM_velp_old, AAM_velp, &
	138	AAM_vel_plus_source(AAM_dissip), AAM_vel_plus_source(AAM_dyn), AAM_vel_plus_source(AAM_phys)
	139	WRITE(,'(A,6E12.4)') 'AAM_vel- : time,old,new,dissip,dyn,phys', dtit, AAM_velm_old, AAM_velm, &
	140	AAM_vel_minus_source(AAM_dissip), AAM_vel_minus_source(AAM_dyn), AAM_vel_minus_source(AAM_phys)
	141	WRITE(*,'(A,6E12.4)') 'AAM_dyn budget mass + vel:', AAM_mass_source(AAM_dyn) + AAM_vel_source(AAM_dyn)
[365]	142	AAM_mass_old=AAM_mass
	143	AAM_vel_old=AAM_vel
[649]	144	AAM_velp_old=AAM_velp
	145	AAM_velm_old=AAM_velm
[365]	146	AAM_mass_source=0.
	147	AAM_vel_source=0.
[649]	148	AAM_vel_minus_source=0.
	149	AAM_vel_plus_source=0.
[365]	150	END IF
[160]	151	END IF
[365]	152
[151]	153	END SUBROUTINE check_conserve
[326]	154
	155	!------------------------ Detailed check (Lauritzen et al., 2014) -----------------------------
	156
[365]	157	SUBROUTINE check_conserve_detailed(it,tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis)
[326]	158	USE pression_mod
	159	USE vorticity_mod
	160	USE caldyn_gcm_mod
	161	USE exner_mod
[365]	162	! USE mpipara, ONLY : is_mpi_root, comm_icosa
	163	USE omp_para, ONLY : is_master
	164	INTEGER, INTENT(IN) :: tag
[326]	165	TYPE(t_field),POINTER :: f_ps(:)
	166	TYPE(t_field),POINTER :: f_dps(:)
	167	TYPE(t_field),POINTER :: f_ue(:)
	168	TYPE(t_field),POINTER :: f_theta_rhodz(:)
	169	TYPE(t_field),POINTER :: f_phis(:)
	170	INTEGER, INTENT(IN) ::it
	171
	172	REAL(rstd),POINTER :: p(:,:),rhodz(:,:)
	173	INTEGER::ind,ierr
[649]	174	REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel, ang_velp, ang_velm
[326]	175
	176	IF(check_type == check_detailed) THEN
	177
	178	CALL transfert_request(f_ue,req_e1_vect)
	179	CALL pression(f_ps,f_p)
	180	DO ind=1,ndomain
	181	IF (.NOT. assigned_domain(ind)) CYCLE
	182	CALL swap_dimensions(ind)
	183	CALL swap_geometry(ind)
	184	p=f_p(ind)
	185	rhodz=f_rhodz(ind)
	186	CALL compute_rhodz(p,rhodz)
	187	END DO
	188	CALL exner(f_ps,f_p,f_pks,f_pk)
[649]	189	CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsv)
[326]	190
[365]	191	IF (is_master) THEN
	192	! XXX=mass,vel
	193	! AAM_XXX = value of total AAM_XXX at last check
	194	! ang_XXX = value of total AAM_XXX now
	195	! AAM_XXX_source(tag) = source of AAA_XXX due to tag=dissip,dyn,phys
	196	!WRITE(*,'(A,3E12.2)') 'AAM_mass', tag, ang_mass, AAM_mass, ang_mass-AAM_mass ! FIXME
	197	!WRITE(*,'(A,3E12.2)') 'AAM_vel ', tag, ang_vel, AAM_vel, ang_vel-AAM_vel ! FIXME
	198	AAM_mass_source(tag) = AAM_mass_source(tag) + ang_mass - AAM_mass
	199	AAM_vel_source(tag) = AAM_vel_source(tag) + ang_vel - AAM_vel
[649]	200	AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_velp - AAM_velp
	201	AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_velm - AAM_velm
	202	!if ((ang_vel - AAM_vel) .gt. 0.) then
	203	! AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_vel - AAM_vel
	204	!else
	205	! AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_vel - AAM_vel
	206	!endif
[365]	207	AAM_mass = ang_mass
	208	AAM_vel = ang_vel
[649]	209	AAM_velp = ang_velp
	210	AAM_velm = ang_velm
[326]	211	END IF
	212	END IF
	213	END SUBROUTINE check_conserve_detailed
	214
[149]	215	!---------------------------------------------------------------------
	216
[326]	217	SUBROUTINE global_sum(loc_sum, glob_sum)
	218	USE mpi_mod
	219	USE mpipara
	220	USE transfert_omp_mod
	221	REAL(rstd), INTENT(IN) :: loc_sum
	222	REAL(rstd), INTENT(OUT) :: glob_sum
	223	REAL(rstd) :: sum
	224
	225	CALL reduce_sum_omp(loc_sum, sum)
	226	IF (using_mpi) THEN
	227	!$OMP MASTER
	228	CALL MPI_REDUCE(sum, glob_sum, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
	229	!$OMP END MASTER
	230	ELSE
	231	glob_sum=loc_sum
[365]	232	ENDIF
[326]	233	END SUBROUTINE global_sum
	234
[160]	235	SUBROUTINE check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
	236	USE mpi_mod
	237	USE mpipara
[186]	238	USE transfert_omp_mod
[295]	239	USE omp_para
[149]	240	IMPLICIT NONE
	241	TYPE(t_field),POINTER :: f_ps(:)
	242	TYPE(t_field),POINTER :: f_dps(:)
	243	REAL(rstd),POINTER :: ps(:),dps(:)
[160]	244	REAL(rstd), INTENT(OUT) :: mtot, rmsdpdt
	245
[149]	246	INTEGER :: ind,i,j,ij
[160]	247	REAL :: mloc, rmsloc
[186]	248	REAL :: mloc_mpi, rmsloc_mpi
[149]	249
[160]	250	mloc=0.0; rmsloc=0.0
[151]	251	DO ind=1,ndomain
[295]	252	IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
[149]	253	CALL swap_dimensions(ind)
	254	CALL swap_geometry(ind)
	255	dps=f_dps(ind)
	256	ps=f_ps(ind)
	257	DO j=jj_begin,jj_end
	258	DO i=ii_begin,ii_end
	259	ij=(j-1)*iim+i
	260	IF (domain(ind)%own(i,j)) THEN
[160]	261	mloc=mloc+ps(ij)*Ai(ij)
	262	rmsloc=rmsloc+dps(ij)*dps(ij)
[151]	263	ENDIF
[149]	264	ENDDO
	265	ENDDO
	266	ENDDO
[172]	267
[326]	268	CALL global_sum(mloc, mtot)
	269	CALL global_sum(rmsloc, rmsdpdt)
[149]	270	END SUBROUTINE check_mass_conserve
	271
	272	!---------------------------------------------------------------------
	273
[326]	274	SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
[649]	275	stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot)
[151]	276	USE pression_mod
	277	USE vorticity_mod
[149]	278	TYPE(t_field), POINTER :: f_ue(:)
	279	TYPE(t_field), POINTER :: f_theta_rhodz(:)
	280	TYPE(t_field), POINTER :: f_phis(:)
[649]	281	REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot
[149]	282	REAL(rstd), POINTER :: ue(:,:)
	283	REAL(rstd), POINTER :: p(:,:)
[387]	284	REAL(rstd), POINTER :: theta_rhodz(:,:,:)
[149]	285	REAL(rstd), POINTER :: pk(:,:)
	286	REAL(rstd), POINTER :: phis(:)
	287	REAL(rstd), POINTER :: rhodz(:,:)
[649]	288	REAL(rstd) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv
[149]	289	INTEGER :: ind
	290
[326]	291	e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0
[649]	292	AAM_velp = 0 ; AAM_velm = 0
[149]	293
[151]	294	DO ind=1,ndomain
[186]	295	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	296	CALL swap_dimensions(ind)
	297	CALL swap_geometry(ind)
	298	ue=f_ue(ind)
	299	p=f_p(ind)
	300	rhodz=f_rhodz(ind)
	301	theta_rhodz=f_theta_rhodz(ind)
	302	pk=f_pk(ind)
	303	phis=f_phis(ind)
[387]	304	CALL compute_energy(ind,ue,p,rhodz,theta_rhodz(:,:,1),pk,phis, &
[649]	305	e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv)
[151]	306	END DO
[295]	307
[326]	308	CALL global_sum(e, etot)
	309	CALL global_sum(s, stot)
	310	CALL global_sum(rmsv, rmsvtot)
	311	CALL global_sum(AAM_mass, AAM_mass_tot)
	312	CALL global_sum(AAM_vel, AAM_vel_tot)
[649]	313	CALL global_sum(AAM_velp, AAM_velp_tot)
	314	CALL global_sum(AAM_velm, AAM_velm_tot)
[326]	315	END SUBROUTINE check_energy
[295]	316
[649]	317	SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv)
[326]	318	USE disvert_mod
	319	USE wind_mod
	320	USE omp_para
[149]	321	INTEGER,INTENT(IN)::ind
	322	REAL(rstd),INTENT(IN) :: u(3iimjjm,llm)
	323	REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
	324	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
	325	REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
	326	REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
	327	REAL(rstd),INTENT(IN) :: phis(iim*jjm)
[649]	328	REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv
[295]	329
[599]	330	REAL(rstd) :: ucenter(iim*jjm,llm,3)
[149]	331	REAL(rstd) :: ulon(iim*jjm,llm)
	332	REAL(rstd) :: ulat(iim*jjm,llm)
[326]	333
	334	REAL(rstd) :: mass_ij, theta_ij, KE_ij, rad,radomeg,lat
	335	INTEGER :: i,j,ij,l
	336
[151]	337	CALL compute_wind_centered(u,ucenter)
	338	CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
	339	rad = radius
	340	radomeg = rad*omega
[326]	341
[295]	342	DO l = ll_begin,ll_end
[326]	343	DO j=jj_begin,jj_end
	344	DO i=ii_begin,ii_end
	345	ij=(j-1)*iim+i
	346	IF (domain(ind)%own(i,j)) THEN
	347	lat = lat_i(ij)
	348	mass_ij = rhodz(ij,l)*Ai(ij)
	349	theta_ij = theta_rhodz(ij,l)/rhodz(ij,l)
	350	KE_ij = 1/(4Ai(ij))(le(ij+u_right)de(ij+u_right)u(ij+u_right,l)**2 + &
	351	le(ij+u_rup)de(ij+u_rup)u(ij+u_rup,l)**2 + &
	352	le(ij+u_lup)de(ij+u_lup)u(ij+u_lup,l)**2 + &
	353	le(ij+u_left)de(ij+u_left)u(ij+u_left,l)**2 + &
	354	le(ij+u_ldown)de(ij+u_ldown)u(ij+u_ldown,l)**2 + &
	355	le(ij+u_rdown)de(ij+u_rdown)u(ij+u_rdown,l)**2 )
	356	rmsv = rmsv + 2mass_ijKE_ij
	357	s = s + Ai(ij)*theta_rhodz(ij,l)
	358	e = e + mass_ij(phis(ij)+theta_ijpk(ij,l)+KE_ij)
	359	AAM_mass = AAM_mass + radcos(lat)mass_ijradomegcos(lat)
	360	AAM_vel = AAM_vel + radcos(lat)mass_ij*ulon(ij,l)
[649]	361	if (ulon(ij,l) > 0.) then
	362	AAM_velp = AAM_velp + radcos(lat)mass_ij*ulon(ij,l)
	363	else
	364	AAM_velm = AAM_velm + radcos(lat)mass_ij*ulon(ij,l)
	365	endif
[326]	366	END IF
	367	END DO
	368	END DO
	369	END DO
	370
[365]	371	END SUBROUTINE compute_energy
[149]	372
	373	!---------------------------------------------------------------------
[365]	374
[295]	375	SUBROUTINE check_PV(ztot)
[365]	376	REAL(rstd),INTENT(OUT) :: ztot
	377	REAL(rstd), POINTER :: vort(:,:)
	378	REAL(rstd), POINTER :: rhodz(:,:)
	379	INTEGER :: ind
	380	REAL(rstd) :: z,z_mpi
	381
	382	z=0
	383	DO ind=1,ndomain
[186]	384	IF (.NOT. assigned_domain(ind)) CYCLE
[151]	385	CALL swap_dimensions(ind)
	386	CALL swap_geometry(ind)
	387	vort=f_vort(ind)
	388	rhodz=f_rhodz(ind)
[295]	389	CALL compute_PV(vort,rhodz,z)
[365]	390	ENDDO
	391	CALL global_sum(z, ztot)
[151]	392	END SUBROUTINE check_PV
[149]	393
[295]	394	SUBROUTINE compute_PV(vort,rhodz,z)
[151]	395	USE disvert_mod
[295]	396	USE omp_para
[149]	397	REAL(rstd),INTENT(IN) :: vort(2iimjjm,llm)
	398	REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
[295]	399	REAL(rstd),INTENT(INOUT) :: z
[149]	400	REAL(rstd)::qv1,qv2
	401	REAL(rstd)::hv1,hv2
	402	INTEGER :: i,j,ij,l,ij2
	403
[151]	404	hv1 = 0.0 ; hv2 = 0.0
[149]	405
[295]	406	DO l = ll_begin,ll_end
[149]	407	DO j=jj_begin+1,jj_end-1
[151]	408	DO i=ii_begin+1,ii_end-1
	409	ij=(j-1)*iim+i
	410	hv1 = Riv2(ij,vup) * rhodz(ij,l) &
	411	+ Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l) &
	412	+ Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
	413
	414	qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
	415
	416	hv2 = Riv2(ij,vdown) * rhodz(ij,l) &
	417	+ Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l) &
	418	+ Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
	419
	420	qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
[149]	421
[295]	422	z = z + (qv1qv1hv1 + qv2qv2hv2)
[151]	423
[149]	424	ENDDO
[151]	425	ENDDO
[295]	426
[151]	427	ENDDO
	428
	429	END SUBROUTINE compute_PV
[149]	430	!---------------------------------------------------------------------
	431
[151]	432	SUBROUTINE compute_rhodz(p,rhodz)
[149]	433	REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
	434	REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
	435	REAL(rstd) ::mass(iim*jjm,llm)
	436	INTEGER :: i,j,ij,l
	437
	438	DO l = 1, llm
	439	DO j=jj_begin,jj_end
	440	DO i=ii_begin,ii_end
	441	ij=(j-1)*iim+i
	442	mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
	443	rhodz(ij,l) = mass(ij,l) / Ai(ij)
	444	ENDDO
	445	ENDDO
	446	ENDDO
	447
[151]	448	END SUBROUTINE compute_rhodz
[149]	449	!====================================================================
	450
	451
	452	END MODULE check_conserve_mod

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

Download in other formats: