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

Last change on this file since 25 was 25, checked in by dubos, 12 years ago

Minor changes :
caldyn_sw.f90, advect_tracer.f90
icosa_mod.f90 : added parameters for NCAR test cases needing global scope
guided_mod.f90 : CALL to guided_ncar now takes tt=it*dt instead of it as input

Significant changes :
timeloop_gcm.f90 : re-activated CALL to advection scheme
disvert_ncar.f90,
etat0_ncar.f90
guided_ncar_mod.f90 : simplification, introduced several getin(...), update due to recent changes in advection test cases (deformational flow, Hadley cell)
run_adv.def : new keys, reorganized for legibility

Tests :
icosa_gcm.exe tested with ncar_adv_shape=const and ncar_adv_wind=solid,deform,hadley.
q1=1 maintained to machine accuracy. Surface pressure slightly oscillates as expected.

FIXME : Tests by Sarvesh with revision 24 show incorrect advection of cosine bell by solid-body rotation. Not fixed.
File size: 11.9 KB
Line 
1MODULE advect_tracer_mod
2  USE icosa
3  PRIVATE
4  INTEGER,PARAMETER::iapp_tracvl= 3
5  REAL(rstd),SAVE :: dt
6
7  TYPE(t_field),POINTER :: f_normal(:)
8  TYPE(t_field),POINTER :: f_tangent(:)
9  TYPE(t_field),POINTER :: f_gradq3d(:)
10
11  PUBLIC init_advect_tracer, advect_tracer
12
13CONTAINS
14
15  SUBROUTINE init_advect_tracer(dt_in)
16    USE advect_mod
17    IMPLICIT NONE
18    REAL(rstd),INTENT(IN) :: dt_in
19    REAL(rstd),POINTER :: tangent(:,:)
20    REAL(rstd),POINTER :: normal(:,:)
21    INTEGER :: ind
22
23    dt=dt_in
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
36  END SUBROUTINE init_advect_tracer
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. 
63    !------------------------------------------------------sarvesh
64    IF ( first ) THEN
65       CALL allocate_field(f_rhodz,field_t,type_real,llm) 
66       CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 
67       CALL allocate_field(f_massflxc,field_u,type_real,llm) 
68       CALL allocate_field(f_massflx,field_u,type_real,llm) 
69       CALL allocate_field(f_uc,field_u,type_real,llm) 
70       first = .FALSE.
71    END IF
72
73    DO ind=1,ndomain
74       CALL swap_dimensions(ind)
75       CALL swap_geometry(ind)
76       rhodz=f_rhodz(ind)
77       massflx=f_massflx(ind)
78       ps=f_ps(ind)
79       u=f_u(ind)
80
81       DO l = 1, llm
82          DO j=jj_begin-1,jj_end+1
83             DO i=ii_begin-1,ii_end+1
84                ij=(j-1)*iim+i
85                rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g 
86             ENDDO
87          ENDDO
88       ENDDO
89
90       DO l = 1, llm
91          DO j=jj_begin-1,jj_end+1
92             DO i=ii_begin-1,ii_end+1
93                ij=(j-1)*iim+i
94                massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
95                massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
96                massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
97             ENDDO
98          ENDDO
99       ENDDO
100    ENDDO
101
102    IF ( iadvtr == 0 ) THEN
103       DO ind=1,ndomain         
104          CALL swap_dimensions(ind)
105          CALL swap_geometry(ind)
106          rhodz=f_rhodz(ind) 
107          rhodzm1 = f_rhodzm1(ind) 
108          massflxc = f_massflxc(ind)  ! accumulated mass fluxes
109          uc = f_uc(ind)              ! time-integrated normal velocity to compute back-trajectory (Miura)
110          rhodzm1 = rhodz
111          massflxc = 0.0 
112          uc = 0.0 
113       END DO
114       CALL transfert_request(f_rhodzm1,req_i1)   !---->
115       CALL transfert_request(f_massflxc,req_e1) !---->
116       CALL transfert_request(f_massflxc,req_e1) !------>
117       CALL transfert_request(f_uc,req_e1) !---->
118       CALL transfert_request(f_uc,req_e1) 
119    END IF
120
121    iadvtr = iadvtr + 1 
122
123    DO ind=1,ndomain
124       CALL swap_dimensions(ind)
125       CALL swap_geometry(ind)
126       massflx  = f_massflx(ind)
127       massflxc = f_massflxc(ind) 
128       uc = f_uc(ind) 
129       u  = f_u(ind) 
130       massflxc = massflxc + massflx 
131       uc = uc + u 
132    END DO
133
134    IF ( iadvtr == iapp_tracvl ) THEN
135       PRINT *, 'Advection scheme'
136       bigt = dt*iapp_tracvl
137       DO ind=1,ndomain
138          CALL swap_dimensions(ind)
139          CALL swap_geometry(ind)
140          uc = f_uc(ind) 
141          uc = uc/real(iapp_tracvl) 
142          massflxc = f_massflx(ind)
143          massflxc = massflxc*dt
144          ! now massflx is time-integrated
145       END DO
146
147       CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt) 
148       iadvtr = 0
149    END IF
150  END SUBROUTINE advect_tracer
151
152  SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt)
153    USE field_mod
154    USE domain_mod
155    USE dimensions
156    USE grid_param
157    USE geometry
158    USE metric
159    USE advect_mod
160    IMPLICIT NONE
161
162    TYPE(t_field),POINTER :: f_q(:)
163    TYPE(t_field),POINTER :: f_u(:)
164    TYPE(t_field),POINTER :: f_rhodz(:) 
165    TYPE(t_field),POINTER :: f_massflx(:)
166
167    TYPE(t_field),POINTER,SAVE :: f_wg(:)
168    TYPE(t_field),POINTER,SAVE :: f_zm(:)
169    TYPE(t_field),POINTER,SAVE :: f_zq(:)
170
171    REAL(rstd)::bigt,dt
172    REAL(rstd),POINTER :: q(:,:,:)
173    REAL(rstd),POINTER :: u(:,:)
174    REAL(rstd),POINTER :: rhodz(:,:) 
175    REAL(rstd),POINTER :: massflx(:,:) 
176    REAL(rstd),POINTER :: wg(:,:)
177    REAL(rstd),POINTER :: zq(:,:,:) 
178    REAL(rstd),POINTER :: zm(:,:) 
179    REAL(rstd),POINTER :: tangent(:,:)
180    REAL(rstd),POINTER :: normal(:,:)
181    REAL(rstd),POINTER :: gradq3d(:,:,:)
182
183    REAL(rstd)::pente_max
184    LOGICAL,SAVE::first = .true. 
185    INTEGER :: i,ij,l,j,ind,k
186    REAL(rstd) :: zzpbar, zzw 
187    REAL::qvmax,qvmin 
188
189    IF ( first ) THEN
190       CALL allocate_field(f_wg,field_t,type_real,llm)
191       CALL allocate_field(f_zm,field_t,type_real,llm)
192       CALL allocate_field(f_zq,field_t,type_real,llm,nqtot) 
193       first = .FALSE.
194    END IF
195
196    DO ind=1,ndomain
197       CALL swap_dimensions(ind)
198       CALL swap_geometry(ind)
199       q=f_q(ind) 
200       rhodz=f_rhodz(ind) 
201       zq=f_zq(ind)
202       zm=f_zm(ind) 
203       zm = rhodz ; zq = q 
204       wg = f_wg(ind)
205       wg = 0.0
206       massflx=f_massflx(ind) 
207       CALL advtrac(massflx,wg) ! compute vertical mass fluxes
208    END DO
209
210    !   CALL transfert_request(f_wg,req_i1)
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.