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

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

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