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

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

Correction for dcmip moist physics
Temperature is ouput instead of virtual temperature.

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