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

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

Transport now working again - tested with dcmip4.1.0
File size: 31.1 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    INTEGER :: i,j,ij,l
305    REAL(rstd) :: ww,uu 
306    REAL(rstd) :: delta
307    REAL(rstd) :: etav,hv, du2
308
309    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
310    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function
311    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:)
312    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
313    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux
314    REAL(rstd),ALLOCATABLE,SAVE :: divm(:,:)  ! mass flux divergence
315    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
316   
317    LOGICAL,SAVE :: first=.TRUE.
318!$OMP THREADPRIVATE(first)
319       
320!$OMP BARRIER     
321!$OMP MASTER 
322!    IF (first) THEN
323    ALLOCATE(theta(iim*jjm,llm))    ! potential temperature
324    ALLOCATE(pk(iim*jjm,llm))       ! Exner function
325    ALLOCATE(pks(iim*jjm))
326    ALLOCATE(alpha(iim*jjm,llm))
327    ALLOCATE(beta(iim*jjm,llm))
328    ALLOCATE(phi(iim*jjm,llm))      ! geopotential
329    ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
330    ALLOCATE(divm(iim*jjm,llm))    ! mass flux convergence
331    ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
332
333!!! Compute theta
334   DO l = 1, llm
335!$OMP DO
336      DO j=jj_begin-1,jj_end+1
337         DO i=ii_begin-1,ii_end+1
338            ij=(j-1)*iim+i
339            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
340         ENDDO
341      ENDDO
342   ENDDO
343
344  DO l = 1, llm
345!!!  Compute mass and theta fluxes
346!$OMP DO
347    DO j=jj_begin-1,jj_end+1
348      DO i=ii_begin-1,ii_end+1
349        ij=(j-1)*iim+i
350        hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
351        hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
352        hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
353
354        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
355        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
356        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
357      ENDDO
358    ENDDO
359
360!!! compute horizontal divergence of fluxes
361!$OMP DO
362    DO j=jj_begin,jj_end
363      DO i=ii_begin,ii_end
364        ij=(j-1)*iim+i 
365        ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
366        divm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  &
367                                ne(ij,rup)*hflux(ij+u_rup,l)       +  & 
368                                ne(ij,lup)*hflux(ij+u_lup,l)       +  & 
369                                ne(ij,left)*hflux(ij+u_left,l)     +  & 
370                                ne(ij,ldown)*hflux(ij+u_ldown,l)   +  & 
371                                ne(ij,rdown)*hflux(ij+u_rdown,l))   
372
373        ! signe ? attention d (rho theta dz)
374        ! dtheta_rhodz =  -div(flux.theta)
375        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
376                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  & 
377                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  & 
378                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  & 
379                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  & 
380                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
381      ENDDO
382    ENDDO
383  ENDDO
384   
385!!! cumulate mass flux divergence from top to bottom
386  DO  l = llm-1, 1, -1
387!$OMP DO
388    DO j=jj_begin,jj_end
389      DO i=ii_begin,ii_end
390        ij=(j-1)*iim+i
391        divm(ij,l) = divm(ij,l) + divm(ij,l+1)
392      ENDDO
393    ENDDO
394  ENDDO
395 
396!!! Compute vertical mass flux
397  DO l = 1,llm-1
398!$OMP DO
399    DO j=jj_begin,jj_end
400      DO i=ii_begin,ii_end
401        ij=(j-1)*iim+i
402        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
403        ! => w>0 for upward transport
404        wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 )
405      ENDDO
406    ENDDO
407  ENDDO
408
409! compute dps, set vertical mass flux at the surface to 0
410!$OMP DO
411  DO j=jj_begin,jj_end
412    DO i=ii_begin,ii_end
413      ij=(j-1)*iim+i
414      wflux(ij,1)  = 0.
415      ! dps/dt = -int(div flux)dz
416      dps(ij)=-divm(ij,1) * g 
417    ENDDO
418  ENDDO
419
420!!! Compute potential vorticity (Coriolis) contribution to du
421
422    SELECT CASE(caldyn_conserv)
423    CASE(energy) ! energy-conserving TRiSK
424
425       DO l=1,llm
426          !$OMP DO
427          DO j=jj_begin,jj_end
428             DO i=ii_begin,ii_end
429                ij=(j-1)*iim+i
430
431                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
432                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
433                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
434                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
435                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
436                     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))+         &
437                     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))+         &
438                     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))+         &
439                     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))+             &
440                     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))
441                du(ij+u_right,l) = .5*uu/de(ij+u_right)
442               
443                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
444                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
445                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
446                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
447                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
448                     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)) +          &
449                     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)) +              &
450                     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)) +              &
451                     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)) +            &
452                     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))
453                du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
454
455               
456                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
457                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
458                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
459                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
460                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
461                     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)) +          &
462                     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)) +        &
463                     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)) +      &
464                     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)) +      &
465                     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))
466                du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
467               
468             ENDDO
469          ENDDO
470       ENDDO
471
472    CASE(enstrophy) ! enstrophy-conserving TRiSK
473 
474       DO l=1,llm
475          !$OMP DO
476          DO j=jj_begin,jj_end
477             DO i=ii_begin,ii_end
478                ij=(j-1)*iim+i
479
480                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
481                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
482                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
483                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
484                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
485                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
486                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
487                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
488                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
489                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
490                du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
491               
492               
493                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
494                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
495                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
496                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
497                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
498                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
499                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
500                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
501                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
502                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
503                du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
504
505                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
506                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
507                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
508                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
509                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
510                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
511                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
512                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
513                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
514                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
515                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
516     
517             ENDDO
518          ENDDO
519       ENDDO
520       
521    CASE DEFAULT
522       STOP
523    END SELECT
524 
525!!! Compute Exner function
526!    PRINT *, 'Computing Exner'
527    CALL compute_exner(ps,p,pks,pk,1) 
528
529!!! Compute geopotential
530
531  ! for first layer
532!$OMP DO
533   DO j=jj_begin-1,jj_end+1
534     DO i=ii_begin-1,ii_end+1
535       ij=(j-1)*iim+i
536       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
537     ENDDO
538   ENDDO
539         
540  ! for other layers
541   DO l = 2, llm
542!$OMP DO
543     DO j=jj_begin-1,jj_end+1
544       DO i=ii_begin-1,ii_end+1
545         ij=(j-1)*iim+i
546         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
547                                       * (  pk(ij,l-1) -  pk(ij,l)    )
548       ENDDO
549     ENDDO
550   ENDDO
551
552   
553!!! Compute bernouilli term = Kinetic Energy + geopotential
554  DO l=1,llm
555!$OMP DO
556    DO j=jj_begin,jj_end
557      DO i=ii_begin,ii_end
558        ij=(j-1)*iim+i
559
560        berni(ij,l) =  phi(ij,l) &                                                         
561                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
562                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
563                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
564                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
565                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
566                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
567       
568      ENDDO
569    ENDDO
570  ENDDO
571 
572 
573!!! gradients of Bernoulli and Exner functions
574  DO l=1,llm
575!$OMP DO
576    DO j=jj_begin,jj_end
577      DO i=ii_begin,ii_end
578        ij=(j-1)*iim+i
579       
580        out_u(ij+u_right,l)=  1/de(ij+u_right) * (                              &
581                          0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
582                             *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
583                        + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
584       
585        du(ij+u_right,l) = du(ij+u_right,l) + out_u(ij+u_right,l)
586       
587        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
588                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
589                             *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
590                        + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
591       
592        du(ij+u_lup,l) = du(ij+u_lup,l) + out_u(ij+u_lup,l)
593               
594        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
595                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
596                             *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
597                        + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
598       
599        du(ij+u_ldown,l) = du(ij+u_ldown,l) + out_u(ij+u_ldown,l)
600
601      ENDDO
602    ENDDO
603  ENDDO
604 
605!!! contributions of vertical advection to du, dtheta
606 
607  DO l=1,llm-1
608!$OMP DO
609    DO j=jj_begin,jj_end
610      DO i=ii_begin,ii_end
611        ij=(j-1)*iim+i
612         ! ww>0 <=> upward transport
613
614        ww            =   0.5 * wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
615        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
616        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
617
618        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_right,l+1))
619        uu  = u(ij+u_right,l+1) - u(ij+u_right,l) 
620        du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
621        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)))
622
623        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_lup,l+1))
624        uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l) 
625        du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
626        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)))
627
628        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_ldown,l+1))
629        uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l) 
630        du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
631        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)))
632
633      ENDDO
634    ENDDO
635  ENDDO     
636   
637!!$OMP BARRIER
638!!$OMP MASTER 
639    DEALLOCATE(theta)   ! potential temperature
640    DEALLOCATE(pk)      ! Exner function
641    DEALLOCATE(pks)
642    DEALLOCATE(alpha)
643    DEALLOCATE(beta)
644    DEALLOCATE(phi)     ! geopotential
645    DEALLOCATE(Ftheta)  ! theta flux   
646    DEALLOCATE(divm)    ! mass flux divergence
647    DEALLOCATE(berni)   ! bernouilli term
648!!$OMP END MASTER
649!!$OMP BARRIER                                                     
650  END SUBROUTINE compute_caldyn
651
652!-------------------------------- Diagnostics ----------------------------
653
654  SUBROUTINE check_mass_conservation(f_ps,f_dps)
655  USE icosa
656  USE mpipara
657  IMPLICIT NONE
658    TYPE(t_field),POINTER :: f_ps(:)
659    TYPE(t_field),POINTER :: f_dps(:)
660    REAL(rstd),POINTER :: ps(:)
661    REAL(rstd),POINTER :: dps(:)
662    REAL(rstd) :: mass_tot,dmass_tot
663    INTEGER :: ind,i,j,ij
664   
665    mass_tot=0
666    dmass_tot=0
667   
668    CALL transfert_request(f_dps,req_i1)
669    CALL transfert_request(f_ps,req_i1)
670
671    DO ind=1,ndomain
672      CALL swap_dimensions(ind)
673      CALL swap_geometry(ind)
674
675      ps=f_ps(ind)
676      dps=f_dps(ind)
677
678      DO j=jj_begin,jj_end
679        DO i=ii_begin,ii_end
680          ij=(j-1)*iim+i
681          IF (domain(ind)%own(i,j)) THEN
682            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
683            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
684          ENDIF
685        ENDDO
686      ENDDO
687   
688    ENDDO
689    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
690
691  END SUBROUTINE check_mass_conservation 
692 
693  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
694       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
695    USE icosa
696    USE vorticity_mod
697    USE theta2theta_rhodz_mod
698    USE pression_mod
699    USE omega_mod
700    USE write_field
701    USE vertical_interp_mod
702    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
703         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
704   
705    REAL(rstd) :: out_pression_lev
706    CHARACTER(LEN=255) :: str_pression
707    CHARACTER(LEN=255) :: physics_type
708   
709    out_pression_level=0
710    CALL getin("out_pression_level",out_pression_level) 
711    WRITE(str_pression,*) INT(out_pression_level/100)
712    str_pression=ADJUSTL(str_pression)
713   
714    CALL writefield("ps",f_ps)
715    CALL writefield("dps",f_dps)
716    CALL writefield("phis",f_phis)
717    CALL vorticity(f_u,f_buf_v)
718    CALL writefield("vort",f_buf_v)
719
720    CALL w_omega(f_ps, f_u, f_buf_i)
721    CALL writefield('omega', f_buf_i)
722    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
723      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
724      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
725    ENDIF
726   
727    ! Temperature
728    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
729   
730    CALL getin('physics',physics_type)
731    IF (TRIM(physics_type)=='dcmip') THEN
732      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
733      CALL writefield("T",f_buf1_i)
734      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
735        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
736        CALL writefield("T"//TRIM(str_pression),f_buf_s)
737      ENDIF
738    ELSE
739      CALL writefield("T",f_buf_i)
740      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
741        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
742        CALL writefield("T"//TRIM(str_pression),f_buf_s)
743      ENDIF
744    ENDIF
745   
746    ! velocity components
747    CALL un2ulonlat(f_u, f_buf_i3, f_buf1_i, f_buf2_i)
748    CALL writefield("ulon",f_buf1_i)
749    CALL writefield("ulat",f_buf2_i)
750
751    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
752      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
753      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
754      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
755      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
756    ENDIF
757   
758    ! geopotential
759    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
760    CALL writefield("p",f_buf_p)
761    CALL writefield("phi",f_buf_i)
762    CALL writefield("theta",f_buf1_i) ! potential temperature
763    CALL writefield("pk",f_buf2_i) ! Exner pressure
764 
765 
766  END SUBROUTINE write_output_fields
767 
768  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
769    USE field_mod
770    USE pression_mod
771    USE exner_mod
772    USE geopotential_mod
773    USE theta2theta_rhodz_mod
774    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
775         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
776    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
777         phi(:,:), phis(:), ps(:), pks(:)
778    INTEGER :: ind
779
780    DO ind=1,ndomain
781       CALL swap_dimensions(ind)
782       CALL swap_geometry(ind)
783       ps = f_ps(ind)
784       p  = f_p(ind)
785       CALL compute_pression(ps,p,0)
786       pk = f_pk(ind)
787       pks = f_pks(ind)
788       CALL compute_exner(ps,p,pks,pk,0)
789       theta_rhodz = f_theta_rhodz(ind)
790       theta = f_theta(ind)
791       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
792       phis = f_phis(ind)
793       phi = f_phi(ind)
794       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
795    END DO
796
797  END SUBROUTINE thetarhodz2geopot
798 
799  SUBROUTINE un2ulonlat(f_u, f_u3d, f_ulon, f_ulat)
800    USE field_mod
801    USE wind_mod
802    TYPE(t_field), POINTER :: f_u(:), & ! IN  : normal velocity components on edges
803         f_u3d(:), f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
804    REAL(rstd),POINTER :: u(:,:), u3d(:,:,:), ulon(:,:), ulat(:,:)
805    INTEGER :: ind
806    DO ind=1,ndomain
807       CALL swap_dimensions(ind)
808       CALL swap_geometry(ind)
809       u=f_u(ind)
810       u3d=f_u3d(ind)
811       CALL compute_wind_centered(u,u3d)
812       ulon=f_ulon(ind)
813       ulat=f_ulat(ind)
814       CALL compute_wind_centered_lonlat_compound(u3d, ulon, ulat)
815    END DO
816  END SUBROUTINE un2ulonlat
817
818  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
819  USE icosa
820  IMPLICIT NONE
821    TYPE(t_field), POINTER :: f_TV(:)
822    TYPE(t_field), POINTER :: f_q(:)
823    TYPE(t_field), POINTER :: f_T(:)
824   
825    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
826    INTEGER :: ind
827   
828    DO ind=1,ndomain
829       CALL swap_dimensions(ind)
830       CALL swap_geometry(ind)
831       Tv=f_Tv(ind)
832       q=f_q(ind)
833       T=f_T(ind)
834       T=Tv/(1+0.608*q(:,:,1))
835    END DO
836   
837  END SUBROUTINE Tv2T
838 
839END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.