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

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

Put time variable : dt, itaumax, write_period, itau_out in the time module

YM
File size: 12.0 KB
Line 
1MODULE advect_tracer_mod
2  USE icosa
3  PRIVATE
4  INTEGER,PARAMETER::iapp_tracvl= 3
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.