source: codes/icosagcm/trunk/src/caldyn_gcm.f90 @ 31

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

Implementation of parallelism
implementation of iterative laplacian for dissipation
File size: 23.8 KB
Line 
1MODULE caldyn_gcm_mod
2  USE icosa
3
4  PRIVATE
5  TYPE(t_field),POINTER :: f_out(:)
6  REAL(rstd),POINTER :: out(:,:)
7  TYPE(t_field),POINTER :: f_out_u(:)
8  REAL(rstd),POINTER :: out_u(:,:)
9  TYPE(t_field),POINTER :: f_out_z(:)
10  REAL(rstd),POINTER :: out_z(:,:)
11 
12  INTEGER :: itau_out
13
14  PUBLIC init_caldyn, caldyn
15
16CONTAINS
17 
18  SUBROUTINE init_caldyn(dt)
19  USE icosa
20  IMPLICIT NONE
21    REAL(rstd),INTENT(IN) :: dt
22    INTEGER :: write_period
23
24    CALL allocate_caldyn
25    CALL getin('write_period',write_period)
26   
27    itau_out=INT(write_period/dt)
28   
29    CALL allocate_caldyn
30 
31  END SUBROUTINE init_caldyn
32 
33  SUBROUTINE allocate_caldyn
34  USE icosa
35  IMPLICIT NONE
36
37    CALL allocate_field(f_out,field_t,type_real,llm) 
38    CALL allocate_field(f_out_u,field_u,type_real,llm) 
39    CALL allocate_field(f_out_z,field_z,type_real,llm) 
40   
41  END SUBROUTINE allocate_caldyn
42 
43  SUBROUTINE swap_caldyn(ind)
44  IMPLICIT NONE
45    INTEGER,INTENT(IN) :: ind
46    out=f_out(ind) 
47    out_u=f_out_u(ind) 
48    out_z=f_out_z(ind) 
49     
50  END SUBROUTINE swap_caldyn
51 
52  SUBROUTINE check_mass_conservation(f_ps,f_dps)
53  USE icosa
54  IMPLICIT NONE
55    TYPE(t_field),POINTER :: f_ps(:)
56    TYPE(t_field),POINTER :: f_dps(:)
57    REAL(rstd),POINTER :: ps(:)
58    REAL(rstd),POINTER :: dps(:)
59    REAL(rstd) :: mass_tot,dmass_tot
60    INTEGER :: ind,i,j,ij
61   
62    mass_tot=0
63    dmass_tot=0
64   
65    CALL transfert_request(f_dps,req_i1)
66    CALL transfert_request(f_ps,req_i1)
67
68    DO ind=1,ndomain
69      CALL swap_dimensions(ind)
70      CALL swap_geometry(ind)
71
72      ps=f_ps(ind)
73      dps=f_dps(ind)
74
75      DO j=jj_begin,jj_end
76        DO i=ii_begin,ii_end
77          ij=(j-1)*iim+i
78          IF (domain(ind)%own(i,j)) THEN
79            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
80            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
81          ENDIF
82        ENDDO
83      ENDDO
84   
85    ENDDO
86    PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
87
88  END SUBROUTINE check_mass_conservation
89 
90 
91 
92  SUBROUTINE caldyn(it,f_phis, f_ps, f_theta_rhodz, f_u, f_dps, f_dtheta_rhodz, f_du)
93  USE icosa
94  USE vorticity_mod
95  USE kinetic_mod
96  USE theta2theta_rhodz_mod
97  IMPLICIT NONE
98  INTEGER,INTENT(IN)    :: it
99  TYPE(t_field),POINTER :: f_phis(:)
100  TYPE(t_field),POINTER :: f_ps(:)
101  TYPE(t_field),POINTER :: f_theta_rhodz(:)
102  TYPE(t_field),POINTER :: f_u(:)
103  TYPE(t_field),POINTER :: f_dps(:)
104  TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
105  TYPE(t_field),POINTER :: f_du(:)
106
107  REAL(rstd),POINTER :: phis(:)
108  REAL(rstd),POINTER :: ps(:)
109  REAL(rstd),POINTER :: theta_rhodz(:,:)
110  REAL(rstd),POINTER :: u(:,:)
111  REAL(rstd),POINTER :: dps(:)
112  REAL(rstd),POINTER :: dtheta_rhodz(:,:)
113  REAL(rstd),POINTER :: du(:,:)
114  INTEGER :: ind,ij
115
116 
117    CALL transfert_request(f_phis,req_i1) 
118    CALL transfert_request(f_ps,req_i1) 
119    CALL transfert_request(f_theta_rhodz,req_i1) 
120    CALL transfert_request(f_u,req_e1)
121!    CALL transfert_request(f_u,req_e1)
122   
123
124!    CALL vorticity(f_u,f_out_z)
125   
126    DO ind=1,ndomain
127      CALL swap_dimensions(ind)
128      CALL swap_geometry(ind)
129      CALL swap_caldyn(ind)
130     
131      phis=f_phis(ind)
132      ps=f_ps(ind)
133      theta_rhodz=f_theta_rhodz(ind)
134      u=f_u(ind)
135      dps=f_dps(ind)
136      dtheta_rhodz=f_dtheta_rhodz(ind)
137      du=f_du(ind)
138!      ij=(jj_end-1-1)*iim+ii_begin
139!      PRINT *,"--> ind=",ind,ij
140!      PRINT *,u(ij+u_right,1)
141!      PRINT *,u(ij+u_rup,1)
142!      PRINT *,u(ij+u_lup,1)
143!      PRINT *,u(ij+u_left,1)
144!      PRINT *,u(ij+u_ldown,1)
145!      PRINT *,u(ij+u_rdown,1)
146
147!      ij=(jj_end-1-1)*iim+ii_end
148!      PRINT *,"--> ind=",ind,ij
149!      PRINT *,u(ij+u_right,1)
150!      PRINT *,u(ij+u_rup,1)
151!      PRINT *,u(ij+u_lup,1)
152!      PRINT *,u(ij+u_left,1)
153!      PRINT *,u(ij+u_ldown,1)
154!      PRINT *,u(ij+u_rdown,1)     
155!$OMP PARALLEL DEFAULT(SHARED)
156      CALL compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
157!$OMP END PARALLEL
158    ENDDO
159
160    CALL transfert_request(f_out_u,req_e1)
161!    CALL transfert_request(f_out_u,req_e1)
162
163!    CALL vorticity(f_u,f_out_z)
164
165    IF (mod(it,itau_out)==0 ) THEN
166      CALL writefield("ps",f_ps)
167      CALL writefield("dps",f_dps)
168!      CALL writefield("theta_rhodz",f_theta_rhodz)
169!    CALL kinetic(f_u,f_out)
170!      CALL writefield("Ki",f_out)
171!      CALL writefield("dtheta_rhodz",f_dtheta_rhodz)
172      CALL vorticity(f_u,f_out_z)
173      CALL writefield("vort",f_out_z)
174!      CALL writefield("theta",f_out)
175      CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_out) ;
176      CALL writefield("T",f_out)
177   
178!      CALL writefield("out",f_out)
179!      DO ind=1,ndomain
180!        CALL writefield("Ki",f_out,ind)
181!        CALL writefield("vort",f_out_z,ind)
182!        CALL writefield("dps",f_dps,ind)
183!      ENDDO
184
185    ENDIF
186!    CALL check_mass_conservation(f_ps,f_dps)
187
188  END SUBROUTINE caldyn
189 
190 
191  SUBROUTINE compute_caldyn(phis, ps, theta_rhodz, u, dps, dtheta_rhodz, du)
192  USE icosa
193  USE disvert_mod
194  IMPLICIT NONE
195    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
196    REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm)
197    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
198    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
199    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
200    REAL(rstd),INTENT(OUT):: dtheta_rhodz(iim*jjm,llm)
201    REAL(rstd),INTENT(OUT):: dps(iim*jjm)
202   
203    INTEGER :: i,j,ij,l
204    REAL(rstd) :: ww,uu 
205    REAL(rstd) :: delta
206    REAL(rstd) :: etav,hv   
207
208!   REAL(rstd) :: theta(iim*jjm,llm)  ! potential temperature
209!   REAL(rstd) :: p(iim*jjm,llm+1)  ! pression
210!   REAL(rstd) :: pk(iim*jjm,llm)   ! Exner function
211!   REAL(rstd) :: pks(iim*jjm)
212!! Intermediate variable to compute exner function
213!   REAL(rstd) :: alpha(iim*jjm,llm)
214!   REAL(rstd) :: beta(iim*jjm,llm)
215!!   
216!   REAL(rstd) :: phi(iim*jjm,llm)    ! geopotential
217!   REAL(rstd) :: mass(iim*jjm,llm)   ! mass   
218!   REAL(rstd) :: rhodz(iim*jjm,llm)   ! mass density   
219!   REAL(rstd) :: Fe(3*iim*jjm,llm)   ! mass flux   
220!   REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux   
221!   REAL(rstd) :: convm(iim*jjm,llm)    ! mass flux convergence
222!   REAL(rstd) :: w(iim*jjm,llm)       ! vertical velocity     
223!   REAL(rstd) :: qv(2*iim*jjm,llm)       ! potential velocity 
224!   REAL(rstd) :: berni(iim*jjm,llm)   ! bernouilli term
225   
226    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
227    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)  ! pression
228    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:)   ! Exner function
229    REAL(rstd),ALLOCATABLE,SAVE :: pks(:)
230 ! Intermediate variable to compute exner function
231    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:)
232    REAL(rstd),ALLOCATABLE,SAVE :: beta(:,:)
233 !   
234    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
235    REAL(rstd),ALLOCATABLE,SAVE :: mass(:,:)   ! mass   
236    REAL(rstd),ALLOCATABLE,SAVE :: rhodz(:,:)   ! mass density   
237    REAL(rstd),ALLOCATABLE,SAVE :: Fe(:,:)   ! mass flux   
238    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux   
239    REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)    ! mass flux convergence
240    REAL(rstd),ALLOCATABLE,SAVE :: w(:,:)       ! vertical velocity     
241    REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)       ! potential velocity 
242    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)   ! bernouilli term
243   
244    LOGICAL,SAVE :: first=.TRUE.
245!$OMP THREADPRIVATE(first)
246       
247!$OMP BARRIER     
248!$OMP MASTER 
249    IF (first) THEN
250      ALLOCATE(theta(iim*jjm,llm))  ! potential temperature
251      ALLOCATE(p(iim*jjm,llm+1))  ! pression
252      ALLOCATE(pk(iim*jjm,llm))   ! Exner function
253      ALLOCATE(pks(iim*jjm))
254      ALLOCATE(alpha(iim*jjm,llm))
255      ALLOCATE(beta(iim*jjm,llm))
256      ALLOCATE(phi(iim*jjm,llm))    ! geopotential
257      ALLOCATE(mass(iim*jjm,llm))   ! mass   
258      ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
259      ALLOCATE(Fe(3*iim*jjm,llm))   ! mass flux   
260      ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
261      ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence
262      ALLOCATE(w(iim*jjm,llm))       ! vertical velocity     
263      ALLOCATE(qv(2*iim*jjm,llm))       ! potential velocity 
264      ALLOCATE(berni(iim*jjm,llm))   ! bernouilli term
265      first=.FALSE.
266    ENDIF
267!$OMP END MASTER
268!$OMP BARRIER       
269 !    du(:,:)=0
270!    theta=1e10
271!    p=1e10
272!    pk=1e10
273!    pks=1e10
274!    alpha=1e10
275!    beta=1e10
276!    phi=1e10
277!    mass=1e10
278!    rhodz=1e10
279!    Fe=1e10
280!    Ftheta=1e10
281!    convm=1e10
282!    w=1e10
283!    qv=1e10
284!    berni=1e10
285   
286 !!! Compute pression
287 
288    DO    l    = 1, llm+1
289!$OMP DO
290      DO j=jj_begin-1,jj_end+1
291        DO i=ii_begin-1,ii_end+1
292          ij=(j-1)*iim+i
293          p(ij,l) = ap(l) + bp(l) * ps(ij)
294        ENDDO
295      ENDDO
296    ENDDO
297     
298
299
300 
301 !!! Compute Exnher function
302 
303 !! Compute Alpha and Beta
304 
305 ! for llm layer
306!$OMP DO
307   DO j=jj_begin-1,jj_end+1
308     DO i=ii_begin-1,ii_end+1
309       ij=(j-1)*iim+i
310       alpha(ij,llm) = 0.
311       beta (ij,llm) = 1./ (1+ 2*kappa)
312     ENDDO
313   ENDDO
314   
315 ! for other layer   
316   DO l = llm-1 , 2 , -1
317!$OMP DO
318     DO j=jj_begin-1,jj_end+1
319       DO i=ii_begin-1,ii_end+1
320          ij=(j-1)*iim+i
321          delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
322          alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
323          beta (ij,l)  =   p(ij,l  ) / delta   
324       ENDDO
325     ENDDO
326   ENDDO
327   
328 !! Compute pk
329 
330 ! for first layer
331!$OMP DO
332   DO j=jj_begin-1,jj_end+1
333     DO i=ii_begin-1,ii_end+1
334        ij=(j-1)*iim+i
335        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
336        pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /               &
337                   (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) )
338     ENDDO
339   ENDDO
340   
341 ! for other layers
342     
343   DO l = 2, llm
344!$OMP DO
345     DO j=jj_begin-1,jj_end+1
346       DO i=ii_begin-1,ii_end+1
347         ij=(j-1)*iim+i
348         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
349       ENDDO
350     ENDDO
351   ENDDO
352   
353 
354!!! Compute mass
355   DO l = 1, llm
356!$OMP DO
357     DO j=jj_begin-1,jj_end+1
358       DO i=ii_begin-1,ii_end+1
359         ij=(j-1)*iim+i
360         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
361         rhodz(ij,l) = mass(ij,l) / Ai(ij)
362       ENDDO
363     ENDDO
364   ENDDO
365
366!! compute theta
367    DO    l    = 1, llm
368!$OMP DO
369      DO j=jj_begin-1,jj_end+1
370        DO i=ii_begin-1,ii_end+1
371          ij=(j-1)*iim+i
372          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
373        ENDDO
374      ENDDO
375    ENDDO
376
377
378!!! Compute geopotential
379
380  ! for first layer
381!$OMP DO
382   DO j=jj_begin-1,jj_end+1
383     DO i=ii_begin-1,ii_end+1
384       ij=(j-1)*iim+i
385       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
386     ENDDO
387   ENDDO
388         
389  ! for other layers
390   DO l = 2, llm
391!$OMP DO
392     DO j=jj_begin-1,jj_end+1
393       DO i=ii_begin-1,ii_end+1
394         ij=(j-1)*iim+i
395         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
396                                       * (  pk(ij,l-1) -  pk(ij,l)    )
397       ENDDO
398     ENDDO
399   ENDDO
400
401   
402!!!  Compute mass flux
403!! question ￜ thomas : meilleure pondￜration de la masse sur les liens ?
404
405  DO l = 1, llm
406!$OMP DO
407    DO j=jj_begin-1,jj_end+1
408      DO i=ii_begin-1,ii_end+1
409        ij=(j-1)*iim+i
410        Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
411        Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
412        Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
413      ENDDO
414    ENDDO
415  ENDDO
416
417!!! fisrt composante dtheta
418
419 ! Flux on the edge
420  DO l = 1, llm
421!$OMP DO
422    DO j=jj_begin-1,jj_end+1
423      DO i=ii_begin-1,ii_end+1
424        ij=(j-1)*iim+i 
425        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*Fe(ij+u_right,l)
426        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*Fe(ij+u_lup,l)
427        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*Fe(ij+u_ldown,l)
428      ENDDO
429    ENDDO
430  ENDDO
431 
432 
433 ! compute divergence
434  DO l = 1, llm
435!$OMP DO
436    DO j=jj_begin,jj_end
437      DO i=ii_begin,ii_end
438        ij=(j-1)*iim+i 
439! signe ? attention d (rho theta dz)
440        ! dtheta_rhodz =  -div(flux.theta)
441        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
442                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  & 
443                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  & 
444                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  & 
445                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  & 
446                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
447      ENDDO
448    ENDDO
449  ENDDO
450
451
452
453!!! mass flux convergence computation 
454
455 ! horizontal convergence
456  DO l = 1, llm
457!$OMP DO
458    DO j=jj_begin,jj_end
459      DO i=ii_begin,ii_end
460        ij=(j-1)*iim+i
461        ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
462        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
463                                ne(ij,rup)*Fe(ij+u_rup,l)       +  & 
464                                ne(ij,lup)*Fe(ij+u_lup,l)       +  & 
465                                ne(ij,left)*Fe(ij+u_left,l)     +  & 
466                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  & 
467                                ne(ij,rdown)*Fe(ij+u_rdown,l))   
468       ENDDO
469     ENDDO
470   ENDDO
471 
472   
473 ! vertical integration from up to down
474  DO  l = llm-1, 1, -1
475!$OMP DO
476    DO j=jj_begin,jj_end
477      DO i=ii_begin,ii_end
478        ij=(j-1)*iim+i
479        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
480      ENDDO
481    ENDDO
482  ENDDO
483 
484
485  DO l = 1, llm
486!$OMP DO
487   DO j=jj_begin,jj_end
488    DO i=ii_begin,ii_end
489      ij=(j-1)*iim+i
490      out(ij,l)=theta(ij,l)-288
491    ENDDO
492  ENDDO
493 ENDDO
494
495
496!!! Compute dps
497!$OMP DO
498  DO j=jj_begin,jj_end
499    DO i=ii_begin,ii_end
500      ij=(j-1)*iim+i
501      ! dps/dt = -int(div flux)dz
502      dps(ij)=-convm(ij,1) * g 
503    ENDDO
504  ENDDO
505
506
507
508!!! Compute vertical velocity
509  DO l = 1,llm-1
510!$OMP DO
511    DO j=jj_begin,jj_end
512      DO i=ii_begin,ii_end
513        ij=(j-1)*iim+i
514        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
515        ! => w>0 for upward transport
516        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
517      ENDDO
518    ENDDO
519  ENDDO
520
521!$OMP DO
522  DO j=jj_begin,jj_end
523    DO i=ii_begin,ii_end
524      ij=(j-1)*iim+i
525      w(ij,1)  = 0.
526    ENDDO
527  ENDDO
528
529
530!!! Compute potential vorticity
531  DO l = 1,llm
532!$OMP DO
533    DO j=jj_begin-1,jj_end+1
534      DO i=ii_begin-1,ii_end+1
535        ij=(j-1)*iim+i
536       
537          etav= 1./Av(ij+z_up)*(  ne(ij,rup)        * u(ij+u_rup,l)        * de(ij+u_rup)         &
538                                + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
539                                - ne(ij,lup)        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
540
541          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
542              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
543              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
544     
545          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
546         
547          etav = 1./Av(ij+z_down)*(  ne(ij,ldown)         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
548                                   + ne(ij+t_ldown,right) * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
549                                   - ne(ij,rdown)         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
550       
551          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
552              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
553              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
554                       
555          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
556
557        ENDDO
558      ENDDO
559    ENDDO
560
561!!! Compute potential vorticity contribution to du
562  DO l=1,llm
563!$OMP DO
564    DO j=jj_begin,jj_end
565      DO i=ii_begin,ii_end
566        ij=(j-1)*iim+i
567
568       du(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))/de(ij+u_right) *      & 
569                      ( wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+                           &
570                        wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+                           &
571                        wee(ij+u_right,3,1)*Fe(ij+u_left,l)+                          &
572                        wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+                         &
573                        wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+                         & 
574                        wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+                 &
575                        wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+                 &
576                        wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+                 &
577                        wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+                   &
578                        wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l) )
579
580
581       du(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))/de(ij+u_lup) *         & 
582                      ( wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+                        &
583                        wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+                       &
584                        wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+                       &
585                        wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+                       &
586                        wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+                         & 
587                        wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+                 &
588                        wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+                   &
589                        wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+                   &
590                        wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+                  &
591                        wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l) )
592
593
594       du(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))/de(ij+u_ldown) *   & 
595                      ( wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+                      &
596                        wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+                      &
597                        wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+                        &
598                        wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+                        &
599                        wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+                       & 
600                        wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+                &
601                        wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+               &
602                        wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+              &
603                        wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+              &
604                        wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l) )
605
606     
607      ENDDO
608    ENDDO
609  ENDDO
610 
611
612!!! Compute bernouilli term = Kinetic Energy + geopotential
613  DO l=1,llm
614!$OMP DO
615    DO j=jj_begin,jj_end
616      DO i=ii_begin,ii_end
617        ij=(j-1)*iim+i
618
619        berni(ij,l) =  phi(ij,l) &                                                         
620                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
621                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
622                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
623                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
624                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
625                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
626       
627      ENDDO
628    ENDDO
629  ENDDO
630 
631 
632!!! second contribution to du
633  DO l=1,llm
634!$OMP DO
635    DO j=jj_begin,jj_end
636      DO i=ii_begin,ii_end
637        ij=(j-1)*iim+i
638       
639         du(ij+u_right,l)= du(ij+u_right,l)+ 1/de(ij+u_right) * (                              &
640                           0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
641                              *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
642                         + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
643         
644         du(ij+u_lup,l)= du(ij+u_lup,l)+ 1/de(ij+u_lup) * (                                  &
645                           0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
646                              *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
647                         + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
648                         
649         du(ij+u_ldown,l)= du(ij+u_ldown,l)+ 1/de(ij+u_ldown) * (                                &
650                           0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
651                              *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
652                         + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
653      ENDDO
654    ENDDO
655  ENDDO
656 
657!!! second contribution to du
658  DO l=1,llm
659!$OMP DO
660    DO j=jj_begin,jj_end
661      DO i=ii_begin,ii_end
662        ij=(j-1)*iim+i
663       
664        out_u(ij+u_right,l)=  1/de(ij+u_right) * (                              &
665                          0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
666                             *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
667                        + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
668       
669        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
670                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
671                             *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
672                        + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
673                       
674        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
675                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
676                             *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
677                        + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
678      ENDDO
679    ENDDO
680  ENDDO
681!!! contribution due to vertical advection
682 
683 ! Contribution to dtheta
684  DO l=1,llm-1
685!$OMP DO
686    DO j=jj_begin,jj_end
687      DO i=ii_begin,ii_end
688         ! ww>0 <=> upward transport
689        ij=(j-1)*iim+i
690        ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
691        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
692        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
693      ENDDO
694    ENDDO
695  ENDDO       
696
697
698 ! Contribution to du
699  DO l=1,llm-1
700!$OMP DO
701    DO j=jj_begin,jj_end
702      DO i=ii_begin,ii_end
703        ij=(j-1)*iim+i
704        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1))
705        uu  = u(ij+u_right,l+1) - u(ij+u_right,l) 
706        du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
707        du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1)))
708
709        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_lup,l+1))
710        uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l) 
711        du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
712        du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1)))
713
714        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_ldown,l+1))
715        uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l) 
716        du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
717        du(ij+u_ldown, l+1 ) = du(ij+u_ldown,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_ldown,l+1)))
718
719      ENDDO
720    ENDDO
721  ENDDO     
722   
723!!$OMP BARRIER
724!!$OMP MASTER 
725!    DEALLOCATE(theta)  ! potential temperature
726!    DEALLOCATE(p)  ! pression
727!    DEALLOCATE(pk)   ! Exner function
728!    DEALLOCATE(pks)
729!    DEALLOCATE(alpha)
730!    DEALLOCATE(beta)
731!    DEALLOCATE(phi)    ! geopotential
732!    DEALLOCATE(mass)   ! mass   
733!    DEALLOCATE(rhodz)   ! mass density   
734!    DEALLOCATE(Fe)   ! mass flux   
735!    DEALLOCATE(Ftheta) ! theta flux   
736!    DEALLOCATE(convm)    ! mass flux convergence
737!    DEALLOCATE(w)       ! vertical velocity     
738!    DEALLOCATE(qv)       ! potential velocity 
739!    DEALLOCATE(berni)   ! bernouilli term
740!!$OMP END MASTER
741!!$OMP BARRIER                                                     
742  END SUBROUTINE compute_caldyn
743 
744 
745 
746END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.