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

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

Simplify the management of the module.

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