source: codes/icosagcm/trunk/src/caldyn_gcm_opt.f90 @ 150

Last change on this file since 150 was 148, checked in by ymipsl, 11 years ago

Various optimisations

  • dissipation is not called every timestep (similar way than LMDZ)
  • transfert size of halos has been reduced : need to synchronise redondant data between tiles at itau_sync timestep

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