Context Navigation

source: codes/icosagcm/trunk/src/advect_tracer.f90 @ 78

Last change on this file since 78 was 30, checked in by sdubey, 12 years ago
bug fix for advection, with transfer
File size: 12.0 KB

Rev	Line
[17]	1	MODULE advect_tracer_mod
[19]	2	USE icosa
[17]	3	PRIVATE
	4	INTEGER,PARAMETER::iapp_tracvl= 3
	5	REAL(rstd),SAVE :: dt
[22]	6
	7	TYPE(t_field),POINTER :: f_normal(:)
	8	TYPE(t_field),POINTER :: f_tangent(:)
	9	TYPE(t_field),POINTER :: f_gradq3d(:)
	10
[17]	11	PUBLIC init_advect_tracer, advect_tracer
	12
	13	CONTAINS
[22]	14
[17]	15	SUBROUTINE init_advect_tracer(dt_in)
[22]	16	USE advect_mod
	17	IMPLICIT NONE
[17]	18	REAL(rstd),INTENT(IN) :: dt_in
[22]	19	REAL(rstd),POINTER :: tangent(:,:)
	20	REAL(rstd),POINTER :: normal(:,:)
[23]	21	INTEGER :: ind
[22]	22
[17]	23	dt=dt_in
[22]	24	CALL allocate_field(f_normal,field_u,type_real,3)
	25	CALL allocate_field(f_tangent,field_u,type_real,3)
	26	CALL allocate_field(f_gradq3d,field_t,type_real,llm,3)
	27
	28	DO ind=1,ndomain
	29	CALL swap_dimensions(ind)
	30	CALL swap_geometry(ind)
	31	normal=f_normal(ind)
	32	tangent=f_tangent(ind)
	33	CALL init_advect(normal,tangent)
	34	END DO
	35
[17]	36	END SUBROUTINE init_advect_tracer
[22]	37
	38	SUBROUTINE advect_tracer(f_ps,f_u,f_q)
	39	USE icosa
	40	USE advect_mod
	41	USE disvert_mod
	42	IMPLICIT NONE
	43	TYPE(t_field),POINTER :: f_ps(:)
	44	TYPE(t_field),POINTER :: f_u(:)
	45	TYPE(t_field),POINTER :: f_q(:)
	46	REAL(rstd),POINTER :: q(:,:,:)
	47	REAL(rstd),POINTER :: u(:,:)
	48	REAL(rstd),POINTER :: ps(:)
	49	REAL(rstd),POINTER :: massflx(:,:)
	50	REAL(rstd),POINTER :: rhodz(:,:)
	51	TYPE(t_field),POINTER,SAVE :: f_massflxc(:)
	52	TYPE(t_field),POINTER,SAVE :: f_massflx(:)
	53	TYPE(t_field),POINTER,SAVE :: f_uc(:)
	54	TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:)
	55	TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
	56	REAL(rstd),POINTER,SAVE :: massflxc(:,:)
	57	REAL(rstd),POINTER,SAVE :: uc(:,:)
	58	REAL(rstd),POINTER,SAVE :: rhodzm1(:,:)
	59	REAL(rstd):: bigt
	60	INTEGER :: ind,it,i,j,l,ij
	61	INTEGER,SAVE :: iadvtr=0
	62	LOGICAL,SAVE:: first=.TRUE.
[30]	63	!------------------------------------------------------sarvesh
	64	CALL transfert_request(f_ps,req_i1)
	65	CALL transfert_request(f_u,req_e1)
	66	CALL transfert_request(f_u,req_e1)
	67	CALL transfert_request(f_q,req_i1)
	68
[17]	69	IF ( first ) THEN
[22]	70	CALL allocate_field(f_rhodz,field_t,type_real,llm)
	71	CALL allocate_field(f_rhodzm1,field_t,type_real,llm)
	72	CALL allocate_field(f_massflxc,field_u,type_real,llm)
	73	CALL allocate_field(f_massflx,field_u,type_real,llm)
	74	CALL allocate_field(f_uc,field_u,type_real,llm)
	75	first = .FALSE.
	76	END IF
[17]	77
	78	DO ind=1,ndomain
[22]	79	CALL swap_dimensions(ind)
	80	CALL swap_geometry(ind)
	81	rhodz=f_rhodz(ind)
	82	massflx=f_massflx(ind)
	83	ps=f_ps(ind)
	84	u=f_u(ind)
	85
	86	DO l = 1, llm
	87	DO j=jj_begin-1,jj_end+1
	88	DO i=ii_begin-1,ii_end+1
	89	ij=(j-1)*iim+i
	90	rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g
	91	ENDDO
[17]	92	ENDDO
[22]	93	ENDDO
[17]	94
[22]	95	DO l = 1, llm
	96	DO j=jj_begin-1,jj_end+1
	97	DO i=ii_begin-1,ii_end+1
	98	ij=(j-1)*iim+i
	99	massflx(ij+u_right,l)=0.5(rhodz(ij,l)+rhodz(ij+t_right,l))u(ij+u_right,l)*le(ij+u_right)
	100	massflx(ij+u_lup,l)=0.5(rhodz(ij,l)+rhodz(ij+t_lup,l))u(ij+u_lup,l)*le(ij+u_lup)
	101	massflx(ij+u_ldown,l)=0.5(rhodz(ij,l)+rhodz(ij+t_ldown,l))u(ij+u_ldown,l)*le(ij+u_ldown)
	102	ENDDO
[17]	103	ENDDO
[22]	104	ENDDO
[17]	105	ENDDO
	106
[22]	107	IF ( iadvtr == 0 ) THEN
	108	DO ind=1,ndomain
	109	CALL swap_dimensions(ind)
	110	CALL swap_geometry(ind)
	111	rhodz=f_rhodz(ind)
	112	rhodzm1 = f_rhodzm1(ind)
	113	massflxc = f_massflxc(ind) ! accumulated mass fluxes
	114	uc = f_uc(ind) ! time-integrated normal velocity to compute back-trajectory (Miura)
	115	rhodzm1 = rhodz
	116	massflxc = 0.0
	117	uc = 0.0
	118	END DO
	119	CALL transfert_request(f_rhodzm1,req_i1) !---->
	120	CALL transfert_request(f_massflxc,req_e1) !---->
	121	CALL transfert_request(f_massflxc,req_e1) !------>
	122	CALL transfert_request(f_uc,req_e1) !---->
	123	CALL transfert_request(f_uc,req_e1)
	124	END IF
[17]	125
[22]	126	iadvtr = iadvtr + 1
[17]	127
[22]	128	DO ind=1,ndomain
[17]	129	CALL swap_dimensions(ind)
	130	CALL swap_geometry(ind)
[22]	131	massflx = f_massflx(ind)
	132	massflxc = f_massflxc(ind)
	133	uc = f_uc(ind)
	134	u = f_u(ind)
	135	massflxc = massflxc + massflx
	136	uc = uc + u
	137	END DO
	138
	139	IF ( iadvtr == iapp_tracvl ) THEN
[25]	140	PRINT *, 'Advection scheme'
[22]	141	bigt = dt*iapp_tracvl
	142	DO ind=1,ndomain
	143	CALL swap_dimensions(ind)
	144	CALL swap_geometry(ind)
	145	uc = f_uc(ind)
	146	uc = uc/real(iapp_tracvl)
[30]	147	massflxc = f_massflxc(ind)
[22]	148	massflxc = massflxc*dt
	149	! now massflx is time-integrated
[17]	150	END DO
	151
	152	CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt)
[22]	153	iadvtr = 0
	154	END IF
[17]	155	END SUBROUTINE advect_tracer
	156
[22]	157	SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt)
	158	USE field_mod
	159	USE domain_mod
	160	USE dimensions
	161	USE grid_param
	162	USE geometry
	163	USE metric
	164	USE advect_mod
	165	IMPLICIT NONE
[17]	166
[22]	167	TYPE(t_field),POINTER :: f_q(:)
	168	TYPE(t_field),POINTER :: f_u(:)
	169	TYPE(t_field),POINTER :: f_rhodz(:)
	170	TYPE(t_field),POINTER :: f_massflx(:)
[17]	171
[22]	172	TYPE(t_field),POINTER,SAVE :: f_wg(:)
	173	TYPE(t_field),POINTER,SAVE :: f_zm(:)
	174	TYPE(t_field),POINTER,SAVE :: f_zq(:)
[17]	175
[22]	176	REAL(rstd)::bigt,dt
	177	REAL(rstd),POINTER :: q(:,:,:)
	178	REAL(rstd),POINTER :: u(:,:)
	179	REAL(rstd),POINTER :: rhodz(:,:)
	180	REAL(rstd),POINTER :: massflx(:,:)
	181	REAL(rstd),POINTER :: wg(:,:)
	182	REAL(rstd),POINTER :: zq(:,:,:)
	183	REAL(rstd),POINTER :: zm(:,:)
	184	REAL(rstd),POINTER :: tangent(:,:)
	185	REAL(rstd),POINTER :: normal(:,:)
	186	REAL(rstd),POINTER :: gradq3d(:,:,:)
[17]	187
[22]	188	REAL(rstd)::pente_max
	189	LOGICAL,SAVE::first = .true.
	190	INTEGER :: i,ij,l,j,ind,k
	191	REAL(rstd) :: zzpbar, zzw
	192	REAL::qvmax,qvmin
[17]	193
[22]	194	IF ( first ) THEN
[17]	195	CALL allocate_field(f_wg,field_t,type_real,llm)
	196	CALL allocate_field(f_zm,field_t,type_real,llm)
	197	CALL allocate_field(f_zq,field_t,type_real,llm,nqtot)
	198	first = .FALSE.
[22]	199	END IF
[17]	200
[22]	201	DO ind=1,ndomain
[17]	202	CALL swap_dimensions(ind)
	203	CALL swap_geometry(ind)
[22]	204	q=f_q(ind)
	205	rhodz=f_rhodz(ind)
	206	zq=f_zq(ind)
	207	zm=f_zm(ind)
	208	zm = rhodz ; zq = q
	209	wg = f_wg(ind)
	210	wg = 0.0
	211	massflx=f_massflx(ind)
	212	CALL advtrac(massflx,wg) ! compute vertical mass fluxes
	213	END DO
[17]	214
[22]	215	DO ind=1,ndomain
	216	CALL swap_dimensions(ind)
	217	CALL swap_geometry(ind)
	218	zq=f_zq(ind)
	219	zm=f_zm(ind)
	220	wg=f_wg(ind)
	221	wg=wg*0.5
[17]	222	DO k = 1, nqtot
[22]	223	CALL vlz(zq(:,:,k),2.,zm,wg)
[17]	224	END DO
[22]	225	END DO
[17]	226
[22]	227	DO ind=1,ndomain
	228	CALL swap_dimensions(ind)
	229	CALL swap_geometry(ind)
	230	zq=f_zq(ind)
	231	zq = f_zq(ind)
	232	zm = f_zm(ind)
	233	massflx =f_massflx(ind)
	234	u = f_u(ind)
	235	tangent = f_tangent(ind)
	236	normal = f_normal(ind)
	237	gradq3d = f_gradq3d(ind)
[17]	238
[22]	239	DO k = 1,nqtot
	240	CALL compute_gradq3d(zq(:,:,k),gradq3d)
[23]	241	CALL compute_advect_horiz(tangent,normal,zq(:,:,k),gradq3d,zm,u,massflx,bigt)
[22]	242	END DO
	243	END DO
[17]	244
[22]	245	DO ind=1,ndomain
	246	CALL swap_dimensions(ind)
	247	CALL swap_geometry(ind)
	248	q = f_q(ind)
	249	zq = f_zq(ind)
	250	zm = f_zm(ind)
	251	wg = f_wg(ind)
	252	DO k = 1,nqtot
	253	CALL vlz(zq(:,:,k),2.,zm,wg)
	254	END DO
	255	q = zq
	256	END DO
[17]	257
[22]	258	END SUBROUTINE vlsplt
[17]	259
[22]	260	!==============================================================================
	261	SUBROUTINE advtrac(massflx,wgg)
	262	USE domain_mod
	263	USE dimensions
	264	USE grid_param
	265	USE geometry
	266	USE metric
	267	USE disvert_mod
	268	IMPLICIT NONE
	269	REAL(rstd),INTENT(IN) :: massflx(iim3jjm,llm)
	270	REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm)
[17]	271
[22]	272	INTEGER :: i,j,ij,l
	273	REAL(rstd) :: convm(iim*jjm,llm)
[17]	274
[22]	275	DO l = 1, llm
	276	DO j=jj_begin,jj_end
	277	DO i=ii_begin,ii_end
	278	ij=(j-1)*iim+i
	279	! divergence of horizontal flux
	280	convm(ij,l)= 1/(Ai(ij))(ne(ij,right)massflx(ij+u_right,l) + &
	281	ne(ij,rup)*massflx(ij+u_rup,l) + &
	282	ne(ij,lup)*massflx(ij+u_lup,l) + &
	283	ne(ij,left)*massflx(ij+u_left,l) + &
	284	ne(ij,ldown)*massflx(ij+u_ldown,l) + &
	285	ne(ij,rdown)*massflx(ij+u_rdown,l))
	286	ENDDO
	287	ENDDO
	288	ENDDO
[17]	289
[22]	290	! accumulate divergence from top of model
	291	DO l = llm-1, 1, -1
	292	DO j=jj_begin,jj_end
	293	DO i=ii_begin,ii_end
	294	ij=(j-1)*iim+i
	295	convm(ij,l) = convm(ij,l) + convm(ij,l+1)
	296	ENDDO
	297	ENDDO
	298	ENDDO
[17]	299
[22]	300	!!! Compute vertical velocity
	301	DO l = 1,llm-1
	302	DO j=jj_begin,jj_end
	303	DO i=ii_begin,ii_end
	304	ij=(j-1)*iim+i
	305	wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ))
	306	ENDDO
	307	ENDDO
	308	ENDDO
[17]	309
[22]	310	DO j=jj_begin,jj_end
	311	DO i=ii_begin,ii_end
	312	ij=(j-1)*iim+i
	313	wgg(ij,1) = 0.
[17]	314	ENDDO
[22]	315	ENDDO
	316	END SUBROUTINE advtrac
[17]	317
[22]	318	SUBROUTINE vlz(q,pente_max,masse,wgg)
	319	!c
	320	!c Auteurs: P.Le Van, F.Hourdin, F.Forget
	321	!c
	322	!c ********************************************************************
	323	!c Shema d'advection " pseudo amont " .
	324	!c ********************************************************************
	325	USE icosa
	326	IMPLICIT NONE
	327	!c
	328	!c Arguments:
	329	!c ----------
	330	REAL masse(iim*jjm,llm),pente_max
	331	REAL q(iim*jjm,llm)
	332	REAL wgg(iimjjm,llm),w(iimjjm,llm+1)
	333	REAL dq(iim*jjm,llm)
	334	INTEGER i,ij,l,j,ii
	335	!c
	336	REAL wq(iim*jjm,llm+1),newmasse
	337
	338	REAL dzq(iimjjm,llm),dzqw(iimjjm,llm),adzqw(iim*jjm,llm),dzqmax
	339	REAL sigw
	340
	341	REAL SSUM
	342
	343
	344	w(:,1:llm) = -wgg(:,:) ! w>0 for downward transport
	345	w(:,llm+1) = 0.0
	346
	347	!c On oriente tout dans le sens de la pression c'est a dire dans le
	348	!c sens de W
	349
	350	DO l=2,llm
	351	DO j=jj_begin,jj_end
	352	DO i=ii_begin,ii_end
	353	ij=(j-1)*iim+i
	354	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
	355	adzqw(ij,l)=abs(dzqw(ij,l))
	356	ENDDO
	357	ENDDO
	358	ENDDO
	359
	360	DO l=2,llm-1
	361	DO j=jj_begin,jj_end
	362	DO i=ii_begin,ii_end
	363	ij=(j-1)*iim+i
	364	IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
[17]	365	dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
[22]	366	ELSE
[17]	367	dzq(ij,l)=0.
[22]	368	ENDIF
	369	dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
	370	dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
	371	ENDDO
	372	ENDDO
	373	ENDDO
[17]	374
[22]	375	DO l=2,llm-1
	376	DO j=jj_begin,jj_end
	377	DO i=ii_begin,ii_end
	378	ij=(j-1)*iim+i
	379	dzq(ij,1)=0.
	380	dzq(ij,llm)=0.
	381	ENDDO
	382	ENDDO
	383	ENDDO
	384	!c ---------------------------------------------------------------
	385	!c .... calcul des termes d'advection verticale .......
	386	!c ---------------------------------------------------------------
[17]	387
[22]	388	!c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
[17]	389
[22]	390	DO l = 1,llm-1
	391	DO j=jj_begin,jj_end
	392	DO i=ii_begin,ii_end
[17]	393	ij=(j-1)*iim+i
	394	IF(w(ij,l+1).gt.0.) THEN
[22]	395	! upwind only if downward transport
	396	sigw=w(ij,l+1)/masse(ij,l+1)
	397	wq(ij,l+1)=w(ij,l+1)(q(ij,l+1)+0.5(1.-sigw)*dzq(ij,l+1))
[17]	398	ELSE
[22]	399	! upwind only if upward transport
	400	sigw=w(ij,l+1)/masse(ij,l)
	401	wq(ij,l+1)=w(ij,l+1)(q(ij,l)-0.5(1.+sigw)*dzq(ij,l))
	402	ENDIF
	403	ENDDO
	404	ENDDO
	405	END DO
[17]	406
[22]	407	DO j=jj_begin,jj_end
	408	DO i=ii_begin,ii_end
	409	ij=(j-1)*iim+i
	410	wq(ij,llm+1)=0.
	411	wq(ij,1)=0.
	412	ENDDO
	413	END DO
[17]	414
[22]	415	DO l=1,llm
[17]	416	DO j=jj_begin,jj_end
[22]	417	DO i=ii_begin,ii_end
[17]	418	ij=(j-1)*iim+i
[22]	419	! masse -= dw/dz but w>0 <=> downward
	420	newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l))
	421	! dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>>
	422	q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse
	423	! masse(ij,l)=newmasse
	424	ENDDO
	425	ENDDO
	426	END DO
	427	RETURN
	428	END SUBROUTINE vlz
[17]	429
	430	END MODULE advect_tracer_mod

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

Download in other formats: