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

Last change on this file since 132 was 132, checked in by dubos, 11 years ago

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