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

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

Simplify the management of the module.

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