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

Last change on this file since 110 was 110, checked in by ymipsl, 12 years ago

Correction for dcmip moist physics
Temperature is ouput instead of virtual temperature.

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