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

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

Implementation of mixte parallelism MPI/OpenMP into src directory

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