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

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

correction for compiling with gfortran (line too long)
improvement for splitting domain
Call twice transfert request for field u is no longer necessary

YM
File size: 23.6 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(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
441                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  & 
442                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  & 
443                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  & 
444                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  & 
445                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
446      ENDDO
447    ENDDO
448  ENDDO
449
450
451
452!!! mass flux convergence computation 
453
454 ! horizontal convergence
455  DO l = 1, llm
456!$OMP DO
457    DO j=jj_begin,jj_end
458      DO i=ii_begin,ii_end
459        ij=(j-1)*iim+i
460!signe ?
461        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
462                                ne(ij,rup)*Fe(ij+u_rup,l)       +  & 
463                                ne(ij,lup)*Fe(ij+u_lup,l)       +  & 
464                                ne(ij,left)*Fe(ij+u_left,l)     +  & 
465                                ne(ij,ldown)*Fe(ij+u_ldown,l)   +  & 
466                                ne(ij,rdown)*Fe(ij+u_rdown,l))   
467       ENDDO
468     ENDDO
469   ENDDO
470 
471   
472 ! vertical integration from up to down
473  DO  l = llm-1, 1, -1
474!$OMP DO
475    DO j=jj_begin,jj_end
476      DO i=ii_begin,ii_end
477        ij=(j-1)*iim+i
478        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
479      ENDDO
480    ENDDO
481  ENDDO
482 
483
484  DO l = 1, llm
485!$OMP DO
486   DO j=jj_begin,jj_end
487    DO i=ii_begin,ii_end
488      ij=(j-1)*iim+i
489      out(ij,l)=theta(ij,l)-288
490    ENDDO
491  ENDDO
492 ENDDO
493
494
495!!! Compute dps
496!$OMP DO
497  DO j=jj_begin,jj_end
498    DO i=ii_begin,ii_end
499      ij=(j-1)*iim+i
500      dps(ij)=-convm(ij,1) * g 
501    ENDDO
502  ENDDO
503
504
505
506!!! Compute vertical velocity
507  DO l = 1,llm-1
508!$OMP DO
509    DO j=jj_begin,jj_end
510      DO i=ii_begin,ii_end
511        ij=(j-1)*iim+i
512        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
513      ENDDO
514    ENDDO
515  ENDDO
516
517!$OMP DO
518  DO j=jj_begin,jj_end
519    DO i=ii_begin,ii_end
520      ij=(j-1)*iim+i
521      w(ij,1)  = 0.
522    ENDDO
523  ENDDO
524
525
526!!! Compute potential vorticity
527  DO l = 1,llm
528!$OMP DO
529    DO j=jj_begin-1,jj_end+1
530      DO i=ii_begin-1,ii_end+1
531        ij=(j-1)*iim+i
532       
533          etav= 1./Av(ij+z_up)*(  ne(ij,rup)        * u(ij+u_rup,l)        * de(ij+u_rup)         &
534                                + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
535                                - ne(ij,lup)        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
536
537          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
538              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
539              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
540     
541          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
542         
543          etav = 1./Av(ij+z_down)*(  ne(ij,ldown)         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
544                                   + ne(ij+t_ldown,right) * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
545                                   - ne(ij,rdown)         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
546       
547          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
548              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
549              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
550                       
551          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
552
553        ENDDO
554      ENDDO
555    ENDDO
556
557!!! Compute potential vorticity contribution to du
558  DO l=1,llm
559!$OMP DO
560    DO j=jj_begin,jj_end
561      DO i=ii_begin,ii_end
562        ij=(j-1)*iim+i
563
564       du(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))/de(ij+u_right) *      & 
565                      ( wee(ij+u_right,1,1)*Fe(ij+u_rup,l)+                           &
566                        wee(ij+u_right,2,1)*Fe(ij+u_lup,l)+                           &
567                        wee(ij+u_right,3,1)*Fe(ij+u_left,l)+                          &
568                        wee(ij+u_right,4,1)*Fe(ij+u_ldown,l)+                         &
569                        wee(ij+u_right,5,1)*Fe(ij+u_rdown,l)+                         & 
570                        wee(ij+u_right,1,2)*Fe(ij+t_right+u_ldown,l)+                 &
571                        wee(ij+u_right,2,2)*Fe(ij+t_right+u_rdown,l)+                 &
572                        wee(ij+u_right,3,2)*Fe(ij+t_right+u_right,l)+                 &
573                        wee(ij+u_right,4,2)*Fe(ij+t_right+u_rup,l)+                   &
574                        wee(ij+u_right,5,2)*Fe(ij+t_right+u_lup,l) )
575
576
577       du(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))/de(ij+u_lup) *         & 
578                      ( wee(ij+u_lup,1,1)*Fe(ij+u_left,l)+                        &
579                        wee(ij+u_lup,2,1)*Fe(ij+u_ldown,l)+                       &
580                        wee(ij+u_lup,3,1)*Fe(ij+u_rdown,l)+                       &
581                        wee(ij+u_lup,4,1)*Fe(ij+u_right,l)+                       &
582                        wee(ij+u_lup,5,1)*Fe(ij+u_rup,l)+                         & 
583                        wee(ij+u_lup,1,2)*Fe(ij+t_lup+u_right,l)+                 &
584                        wee(ij+u_lup,2,2)*Fe(ij+t_lup+u_rup,l)+                   &
585                        wee(ij+u_lup,3,2)*Fe(ij+t_lup+u_lup,l)+                   &
586                        wee(ij+u_lup,4,2)*Fe(ij+t_lup+u_left,l)+                  &
587                        wee(ij+u_lup,5,2)*Fe(ij+t_lup+u_ldown,l) )
588
589
590       du(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))/de(ij+u_ldown) *   & 
591                      ( wee(ij+u_ldown,1,1)*Fe(ij+u_rdown,l)+                      &
592                        wee(ij+u_ldown,2,1)*Fe(ij+u_right,l)+                      &
593                        wee(ij+u_ldown,3,1)*Fe(ij+u_rup,l)+                        &
594                        wee(ij+u_ldown,4,1)*Fe(ij+u_lup,l)+                        &
595                        wee(ij+u_ldown,5,1)*Fe(ij+u_left,l)+                       & 
596                        wee(ij+u_ldown,1,2)*Fe(ij+t_ldown+u_lup,l)+                &
597                        wee(ij+u_ldown,2,2)*Fe(ij+t_ldown+u_left,l)+               &
598                        wee(ij+u_ldown,3,2)*Fe(ij+t_ldown+u_ldown,l)+              &
599                        wee(ij+u_ldown,4,2)*Fe(ij+t_ldown+u_rdown,l)+              &
600                        wee(ij+u_ldown,5,2)*Fe(ij+t_ldown+u_right,l) )
601
602     
603      ENDDO
604    ENDDO
605  ENDDO
606 
607
608!!! Compute bernouilli term = Kinetic Energy + geopotential
609  DO l=1,llm
610!$OMP DO
611    DO j=jj_begin,jj_end
612      DO i=ii_begin,ii_end
613        ij=(j-1)*iim+i
614
615        berni(ij,l) =  phi(ij,l) &                                                         
616                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
617                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
618                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
619                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
620                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
621                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
622       
623      ENDDO
624    ENDDO
625  ENDDO
626 
627 
628!!! second contribution to du
629  DO l=1,llm
630!$OMP DO
631    DO j=jj_begin,jj_end
632      DO i=ii_begin,ii_end
633        ij=(j-1)*iim+i
634       
635         du(ij+u_right,l)= du(ij+u_right,l)+ 1/de(ij+u_right) * (                              &
636                           0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
637                              *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
638                         + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
639         
640         du(ij+u_lup,l)= du(ij+u_lup,l)+ 1/de(ij+u_lup) * (                                  &
641                           0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
642                              *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
643                         + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
644                         
645         du(ij+u_ldown,l)= du(ij+u_ldown,l)+ 1/de(ij+u_ldown) * (                                &
646                           0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
647                              *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
648                         + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
649      ENDDO
650    ENDDO
651  ENDDO
652 
653!!! second contribution to du
654  DO l=1,llm
655!$OMP DO
656    DO j=jj_begin,jj_end
657      DO i=ii_begin,ii_end
658        ij=(j-1)*iim+i
659       
660        out_u(ij+u_right,l)=  1/de(ij+u_right) * (                              &
661                          0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
662                             *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
663                        + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
664       
665        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
666                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
667                             *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
668                        + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
669                       
670        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
671                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
672                             *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
673                        + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
674      ENDDO
675    ENDDO
676  ENDDO
677!!! contribution due to vertical advection
678 
679 ! Contribution to dtheta
680  DO l=1,llm-1
681!$OMP DO
682    DO j=jj_begin,jj_end
683      DO i=ii_begin,ii_end
684        ij=(j-1)*iim+i
685        ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
686        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww 
687        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
688      ENDDO
689    ENDDO
690  ENDDO       
691
692
693 ! Contribution to du
694  DO l=1,llm-1
695!$OMP DO
696    DO j=jj_begin,jj_end
697      DO i=ii_begin,ii_end
698        ij=(j-1)*iim+i
699        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_right,l+1))
700        uu  = u(ij+u_right,l+1) - u(ij+u_right,l) 
701        du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
702        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)))
703
704        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_lup,l+1))
705        uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l) 
706        du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
707        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)))
708
709        ww  = 0.5 * ( w(ij,l+1) + w(ij+t_ldown,l+1))
710        uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l) 
711        du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
712        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)))
713
714      ENDDO
715    ENDDO
716  ENDDO     
717   
718!!$OMP BARRIER
719!!$OMP MASTER 
720!    DEALLOCATE(theta)  ! potential temperature
721!    DEALLOCATE(p)  ! pression
722!    DEALLOCATE(pk)   ! Exner function
723!    DEALLOCATE(pks)
724!    DEALLOCATE(alpha)
725!    DEALLOCATE(beta)
726!    DEALLOCATE(phi)    ! geopotential
727!    DEALLOCATE(mass)   ! mass   
728!    DEALLOCATE(rhodz)   ! mass density   
729!    DEALLOCATE(Fe)   ! mass flux   
730!    DEALLOCATE(Ftheta) ! theta flux   
731!    DEALLOCATE(convm)    ! mass flux convergence
732!    DEALLOCATE(w)       ! vertical velocity     
733!    DEALLOCATE(qv)       ! potential velocity 
734!    DEALLOCATE(berni)   ! bernouilli term
735!!$OMP END MASTER
736!!$OMP BARRIER                                                     
737  END SUBROUTINE compute_caldyn
738 
739 
740 
741END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.