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

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

Some operations must be only done by the mpi master task.

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