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

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

bug fix in caldyn_gcm, didn't work with different size of local domain

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