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

Last change on this file since 132 was 132, checked in by dubos, 11 years ago

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