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

Last change on this file since 78 was 56, checked in by dubos, 12 years ago

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