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

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

minor correction to exclude output at pression level 0

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