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

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

externalize earth constant in run.def and add scalling factor functionnality

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