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

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

caldyn_adv is back
File size: 31.0 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=lmdz
44    CASE('direct')
45       caldyn_exner=direct
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=lmdz
56    CASE('direct')
57       caldyn_hydrostat=direct
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, &
87       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
88    USE icosa
89    USE vorticity_mod
90    USE kinetic_mod
91    USE theta2theta_rhodz_mod
92    USE mpipara
93    IMPLICIT NONE
94    LOGICAL,INTENT(IN)    :: write_out
95    TYPE(t_field),POINTER :: f_phis(:)
96    TYPE(t_field),POINTER :: f_ps(:)
97    TYPE(t_field),POINTER :: f_theta_rhodz(:)
98    TYPE(t_field),POINTER :: f_u(:)
99    TYPE(t_field),POINTER :: f_q(:)
100    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
101    TYPE(t_field),POINTER :: f_dps(:)
102    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
103    TYPE(t_field),POINTER :: f_du(:)
104   
105    REAL(rstd),POINTER :: phis(:), ps(:), dps(:)
106    REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
107    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
108    REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
109    INTEGER :: ind,ij
110   
111    CALL transfert_request(f_phis,req_i1) 
112    CALL transfert_request(f_ps,req_i1) 
113    CALL transfert_request(f_theta_rhodz,req_i1) 
114    CALL transfert_request(f_u,req_e1)
115
116    SELECT CASE(caldyn_conserv)
117    CASE(energy) ! energy-conserving
118       DO ind=1,ndomain
119          CALL swap_dimensions(ind)
120          CALL swap_geometry(ind)
121          ps=f_ps(ind)
122          rhodz=f_rhodz(ind)
123          p=f_p(ind)
124          qu=f_qu(ind)
125          u=f_u(ind)
126          !$OMP PARALLEL DEFAULT(SHARED)
127          CALL compute_pvort(ps, u, p,rhodz,qu)
128          !$OMP END PARALLEL
129       ENDDO
130
131       CALL transfert_request(f_qu,req_e1)
132
133       DO ind=1,ndomain
134          CALL swap_dimensions(ind)
135          CALL swap_geometry(ind)
136          phis=f_phis(ind)
137          hflux=f_hflux(ind)
138          wflux=f_wflux(ind)
139          ps=f_ps(ind)
140          dps=f_dps(ind)
141          theta_rhodz=f_theta_rhodz(ind)
142          dtheta_rhodz=f_dtheta_rhodz(ind)
143          rhodz=f_rhodz(ind)
144          p=f_p(ind)
145          qu=f_qu(ind)
146          u=f_u(ind)
147          du=f_du(ind)
148          out_u=f_out_u(ind) 
149          !$OMP PARALLEL DEFAULT(SHARED)
150          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
151               hflux, wflux, dps, dtheta_rhodz, du)
152          !$OMP END PARALLEL
153       ENDDO       
154       
155    CASE(enstrophy) ! enstrophy-conserving
156       DO ind=1,ndomain
157          CALL swap_dimensions(ind)
158          CALL swap_geometry(ind)
159          phis=f_phis(ind)
160          ps=f_ps(ind)
161          dps=f_dps(ind)
162          hflux=f_hflux(ind)
163          wflux=f_wflux(ind)
164          theta_rhodz=f_theta_rhodz(ind)
165          dtheta_rhodz=f_dtheta_rhodz(ind)
166          rhodz=f_rhodz(ind)
167          p=f_p(ind)
168          qu=f_qu(ind)
169          u=f_u(ind)
170          du=f_du(ind)
171          out_u=f_out_u(ind) 
172          !$OMP PARALLEL DEFAULT(SHARED)
173          CALL compute_pvort(ps, u, p,rhodz,qu)
174          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
175               hflux, wflux, dps, dtheta_rhodz, du)
176          !$OMP END PARALLEL
177       ENDDO
178       
179    CASE DEFAULT
180       STOP
181    END SELECT
182
183    IF (write_out) THEN
184       IF (is_mpi_root) PRINT *,'CALL write_output_fields'
185       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
186            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
187    END IF
188   
189    !    CALL check_mass_conservation(f_ps,f_dps)
190   
191END SUBROUTINE caldyn
192   
193SUBROUTINE compute_pvort(ps, u, p,rhodz,qu)
194  USE icosa
195  USE disvert_mod
196  USE exner_mod
197  IMPLICIT NONE
198  REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
199  REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
200  REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1)
201  REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm)
202  REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
203 
204  INTEGER :: i,j,ij,l
205  REAL(rstd) :: etav,hv
206  REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity 
207 
208  LOGICAL,SAVE :: first=.TRUE.
209  !$OMP THREADPRIVATE(first)
210 
211  !$OMP BARRIER     
212  !$OMP MASTER 
213  !    IF (first) THEN
214  ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity 
215 
216!!! Compute pressure
217  DO    l    = 1, llm+1
218     !$OMP DO
219     DO j=jj_begin-1,jj_end+1
220        DO i=ii_begin-1,ii_end+1
221           ij=(j-1)*iim+i
222           p(ij,l) = ap(l) + bp(l) * ps(ij)
223        ENDDO
224     ENDDO
225  ENDDO
226 
227!!! Compute mass
228  DO l = 1, llm
229     !$OMP DO
230     DO j=jj_begin-1,jj_end+1
231        DO i=ii_begin-1,ii_end+1
232           ij=(j-1)*iim+i
233           rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
234        ENDDO
235     ENDDO
236  ENDDO
237 
238!!! Compute shallow-water potential vorticity
239  DO l = 1,llm
240!$OMP DO
241    DO j=jj_begin-1,jj_end+1
242       DO i=ii_begin-1,ii_end+1
243          ij=(j-1)*iim+i
244         
245          etav= 1./Av(ij+z_up)*(  ne(ij,rup)        * u(ij+u_rup,l)        * de(ij+u_rup)         &
246                                + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
247                                - ne(ij,lup)        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
248
249          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
250              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
251              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
252     
253          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
254         
255          etav = 1./Av(ij+z_down)*(  ne(ij,ldown)         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
256                                   + ne(ij+t_ldown,right) * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
257                                   - ne(ij,rdown)         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
258       
259          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
260              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
261              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
262                       
263          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
264
265        ENDDO
266      ENDDO
267
268      DO j=jj_begin,jj_end
269         DO i=ii_begin,ii_end
270            ij=(j-1)*iim+i
271            qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 
272            qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 
273            qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 
274         END DO
275      END DO
276     
277   ENDDO
278
279!!$OMP BARRIER
280!!$OMP MASTER 
281    DEALLOCATE(qv)       ! potential velocity 
282!!$OMP END MASTER
283!!$OMP BARRIER                                                     
284  END SUBROUTINE compute_pvort
285 
286  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du)
287  USE icosa
288  USE disvert_mod
289  USE exner_mod
290  IMPLICIT NONE
291    REAL(rstd),INTENT(IN)  :: phis(iim*jjm)
292    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
293    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
294    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
295    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1)
296    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
297    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
298
299    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s
300    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
301    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
302    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
303   
304    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
305    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function
306    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:)
307    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
308    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux
309    REAL(rstd),ALLOCATABLE,SAVE :: divm(:,:)  ! mass flux divergence
310    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
311   
312    INTEGER :: i,j,ij,l
313    REAL(rstd) :: ww,uu
314
315    LOGICAL,SAVE :: first=.TRUE.
316!$OMP THREADPRIVATE(first)
317       
318!$OMP BARRIER     
319!$OMP MASTER 
320!    IF (first) THEN
321    ALLOCATE(theta(iim*jjm,llm))    ! potential temperature
322    ALLOCATE(pk(iim*jjm,llm))       ! Exner function
323    ALLOCATE(pks(iim*jjm))
324    ALLOCATE(alpha(iim*jjm,llm))
325    ALLOCATE(beta(iim*jjm,llm))
326    ALLOCATE(phi(iim*jjm,llm))      ! geopotential
327    ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
328    ALLOCATE(divm(iim*jjm,llm))    ! mass flux divvergence
329    ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
330
331!!! Compute theta
332   DO l = 1, llm
333!$OMP DO
334      DO j=jj_begin-1,jj_end+1
335         DO i=ii_begin-1,ii_end+1
336            ij=(j-1)*iim+i
337            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
338         ENDDO
339      ENDDO
340   ENDDO
341
342  DO l = 1, llm
343!!!  Compute mass and theta fluxes
344    DO j=jj_begin-1,jj_end+1
345      DO i=ii_begin-1,ii_end+1
346        ij=(j-1)*iim+i
347        hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
348        hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
349        hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
350
351        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
352        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
353        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
354      ENDDO
355    ENDDO
356
357!!! compute horizontal divergence of fluxes
358    DO j=jj_begin,jj_end
359      DO i=ii_begin,ii_end
360        ij=(j-1)*iim+i 
361        ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
362        divm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  &
363                                ne(ij,rup)*hflux(ij+u_rup,l)       +  & 
364                                ne(ij,lup)*hflux(ij+u_lup,l)       +  & 
365                                ne(ij,left)*hflux(ij+u_left,l)     +  & 
366                                ne(ij,ldown)*hflux(ij+u_ldown,l)   +  & 
367                                ne(ij,rdown)*hflux(ij+u_rdown,l))   
368
369        ! signe ? attention d (rho theta dz)
370        ! dtheta_rhodz =  -div(flux.theta)
371        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
372                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  & 
373                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  & 
374                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  & 
375                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  & 
376                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
377      ENDDO
378    ENDDO
379  ENDDO
380   
381!!! cumulate mass flux divergence from top to bottom
382  DO  l = llm-1, 1, -1
383!$OMP DO
384    DO j=jj_begin,jj_end
385      DO i=ii_begin,ii_end
386        ij=(j-1)*iim+i
387        divm(ij,l) = divm(ij,l) + divm(ij,l+1)
388      ENDDO
389    ENDDO
390  ENDDO
391 
392!!! Compute vertical mass flux
393  DO l = 1,llm-1
394!$OMP DO
395    DO j=jj_begin,jj_end
396      DO i=ii_begin,ii_end
397        ij=(j-1)*iim+i
398        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
399        ! => w>0 for upward transport
400        wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 )
401      ENDDO
402    ENDDO
403  ENDDO
404
405! compute dps, set vertical mass flux at the surface to 0
406!$OMP DO
407  DO j=jj_begin,jj_end
408    DO i=ii_begin,ii_end
409      ij=(j-1)*iim+i
410      wflux(ij,1)  = 0.
411      ! dps/dt = -int(div flux)dz
412      dps(ij)=-divm(ij,1) * g 
413    ENDDO
414  ENDDO
415
416!!! Compute potential vorticity (Coriolis) contribution to du
417
418    SELECT CASE(caldyn_conserv)
419    CASE(energy) ! energy-conserving TRiSK
420
421       DO l=1,llm
422          !$OMP DO
423          DO j=jj_begin,jj_end
424             DO i=ii_begin,ii_end
425                ij=(j-1)*iim+i
426
427                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
428                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
429                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
430                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
431                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
432                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
433                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
434                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
435                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
436                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
437                du(ij+u_right,l) = .5*uu/de(ij+u_right)
438               
439                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
440                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
441                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
442                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
443                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
444                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
445                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
446                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
447                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
448                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
449                du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
450
451               
452                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
453                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
454                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
455                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
456                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
457                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
458                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
459                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
460                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
461                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
462                du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
463               
464             ENDDO
465          ENDDO
466       ENDDO
467
468    CASE(enstrophy) ! enstrophy-conserving TRiSK
469 
470       DO l=1,llm
471          !$OMP DO
472          DO j=jj_begin,jj_end
473             DO i=ii_begin,ii_end
474                ij=(j-1)*iim+i
475
476                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
477                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
478                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
479                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
480                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
481                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
482                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
483                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
484                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
485                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
486                du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
487               
488               
489                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
490                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
491                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
492                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
493                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
494                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
495                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
496                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
497                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
498                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
499                du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
500
501                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
502                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
503                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
504                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
505                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
506                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
507                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
508                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
509                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
510                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
511                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
512     
513             ENDDO
514          ENDDO
515       ENDDO
516       
517    CASE DEFAULT
518       STOP
519    END SELECT
520 
521!!! Compute Exner function
522!    PRINT *, 'Computing Exner'
523    CALL compute_exner(ps,p,pks,pk,1) 
524
525!!! Compute geopotential
526
527  ! for first layer
528!$OMP DO
529   DO j=jj_begin-1,jj_end+1
530     DO i=ii_begin-1,ii_end+1
531       ij=(j-1)*iim+i
532       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
533     ENDDO
534   ENDDO
535         
536  ! for other layers
537   DO l = 2, llm
538!$OMP DO
539     DO j=jj_begin-1,jj_end+1
540       DO i=ii_begin-1,ii_end+1
541         ij=(j-1)*iim+i
542         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
543                                       * (  pk(ij,l-1) -  pk(ij,l)    )
544       ENDDO
545     ENDDO
546   ENDDO
547
548   
549!!! Compute bernouilli term = Kinetic Energy + geopotential
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        berni(ij,l) =  phi(ij,l) &                                                         
557                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
558                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
559                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
560                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
561                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
562                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
563       
564      ENDDO
565    ENDDO
566  ENDDO
567 
568 
569!!! gradients of Bernoulli and Exner functions
570  DO l=1,llm
571!$OMP DO
572    DO j=jj_begin,jj_end
573      DO i=ii_begin,ii_end
574        ij=(j-1)*iim+i
575       
576        out_u(ij+u_right,l)=  1/de(ij+u_right) * (                              &
577                          0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
578                             *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
579                        + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
580       
581        du(ij+u_right,l) = du(ij+u_right,l) + out_u(ij+u_right,l)
582       
583        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
584                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
585                             *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
586                        + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
587       
588        du(ij+u_lup,l) = du(ij+u_lup,l) + out_u(ij+u_lup,l)
589               
590        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
591                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
592                             *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
593                        + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
594       
595        du(ij+u_ldown,l) = du(ij+u_ldown,l) + out_u(ij+u_ldown,l)
596
597      ENDDO
598    ENDDO
599  ENDDO
600 
601!!! contributions of vertical advection to du, dtheta
602 
603  DO l=1,llm-1
604!$OMP DO
605    DO j=jj_begin,jj_end
606      DO i=ii_begin,ii_end
607        ij=(j-1)*iim+i
608         ! ww>0 <=> upward transport
609
610        ww            =   0.5 * wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
611        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
612        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
613
614        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_right,l+1))
615        uu  = u(ij+u_right,l+1) - u(ij+u_right,l) 
616        du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
617        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)))
618
619        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_lup,l+1))
620        uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l) 
621        du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
622        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)))
623
624        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_ldown,l+1))
625        uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l) 
626        du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
627        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)))
628
629      ENDDO
630    ENDDO
631  ENDDO     
632   
633!!$OMP BARRIER
634!!$OMP MASTER 
635    DEALLOCATE(theta)   ! potential temperature
636    DEALLOCATE(pk)      ! Exner function
637    DEALLOCATE(pks)
638    DEALLOCATE(alpha)
639    DEALLOCATE(beta)
640    DEALLOCATE(phi)     ! geopotential
641    DEALLOCATE(Ftheta)  ! theta flux   
642    DEALLOCATE(divm)    ! mass flux divergence
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.