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

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

Merge advection scheme from sarvesh in standard version

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