source: codes/icosagcm/trunk/src/diagnostics/check_conserve.f90

Last change on this file was 953, checked in by adurocher, 5 years ago

trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )
File size: 17.2 KB
Line 
1MODULE check_conserve_mod
2  USE icosa 
3  USE abort_mod
4  IMPLICIT NONE
5
6  PRIVATE
7 
8  INTEGER, PARAMETER, PUBLIC :: check_basic=1, check_detailed=2, &
9       AAM_dissip=1, AAM_dyn=2, AAM_phys=3
10  INTEGER, SAVE :: check_type
11
12  TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
13  TYPE(t_field),POINTER,SAVE :: f_vort(:)
14  TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
15 
16  REAL(rstd),SAVE :: AAM_mass, AAM_mass_old, AAM_vel, AAM_vel_old, &
17       AAM_velp, AAM_velp_old, AAM_velm, AAM_velm_old, &
18       AAM_mass_source(3), AAM_vel_source(3) ! read/written only IF is_master
19  REAL(rstd),SAVE :: AAM_vel_plus_source(3), AAM_vel_minus_source(3)
20  REAL(rstd),SAVE :: mtot0,mqtot0,ztot0,etot0,angtot0,stot0
21  !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0)
22 
23  PUBLIC :: init_check_conserve, check_conserve_detailed, check_conserve
24CONTAINS 
25
26!---------------------------------------------------------------------
27
28  SUBROUTINE init_check_conserve
29    USE getin_mod
30    USE omp_para, ONLY : is_master
31    USE abort_mod
32    CHARACTER(LEN=255) :: check_type_str
33    CALL allocate_field(f_pk,field_t,type_real,llm)
34    CALL allocate_field(f_p,field_t,type_real,llm+1)
35    CALL allocate_field(f_pks,field_t,type_real)
36    CALL allocate_field(f_rhodz,field_t,type_real,llm)
37    CALL allocate_field(f_vort,field_z,type_real,llm)
38   
39    check_type_str='basic'
40    CALL getin("check_conservation",check_type_str)
41    SELECT CASE(TRIM(check_type_str))
42    CASE('basic')
43       check_type=check_basic
44    CASE('detailed')
45       check_type=check_detailed
46    CASE DEFAULT
47       IF(is_master) PRINT*,'Bad selector <', TRIM(check_type_str), &
48            '> for variable check_conservation : options are <basic>, <detailed>'
49       STOP
50    END SELECT
51
52    IF (check_type /= check_basic) THEN
53       CALL abort_acc("check_conservation /= 'basic'")
54    END IF
55  END SUBROUTINE init_check_conserve
56
57!--------------------------------- Basic check --------------------------------
58
59  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,f_q,it) 
60  USE pression_mod
61  USE vorticity_mod
62  USE caldyn_gcm_mod
63  USE exner_mod 
64!  USE mpipara, ONLY : is_mpi_master, comm_icosa
65  USE omp_para, ONLY : is_master
66    TYPE(t_field),POINTER :: f_ps(:)
67    TYPE(t_field),POINTER :: f_dps(:)
68    TYPE(t_field),POINTER :: f_ue(:)
69    TYPE(t_field),POINTER :: f_theta_rhodz(:)
70    TYPE(t_field),POINTER :: f_phis(:)
71    TYPE(t_field),POINTER :: f_q(:)
72    INTEGER, INTENT(IN) :: it
73
74    REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
75    INTEGER :: ind
76    REAL(rstd) :: mtot, mqtot, angtot, rmsdpdt
77    REAL(rstd) :: etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot, ztot
78
79    CALL transfert_request(f_ue,req_e1_vect)
80    CALL pression(f_ps,f_p)
81
82    DO ind=1,ndomain 
83      IF (.NOT. assigned_domain(ind)) CYCLE
84      CALL swap_dimensions(ind)
85      CALL swap_geometry(ind)
86      p=f_p(ind) 
87      rhodz=f_rhodz(ind)
88      CALL compute_rhodz(p,rhodz) 
89    END DO
90
91    CALL vorticity(f_ue,f_vort)
92    CALL check_qmass_conserve(f_q,f_rhodz,mqtot)
93    CALL check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
94    CALL check_PV(ztot) 
95    CALL exner(f_ps,f_p,f_pks,f_pk)
96    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot) 
97
98    IF (is_master) THEN ! is_master = is_omp_master && is_mpi_master
99       ! AAM_mass and AAM_vel must be read/written only IF is_master
100       AAM_mass = ang_mass
101       AAM_vel = ang_vel
102       AAM_velp = ang_velp
103       AAM_velm = ang_velm
104       angtot  = ang_mass + ang_vel
105       IF ( it == itau0  ) THEN
106          ztot0 = ztot
107          mtot0 = mtot
108          mqtot0 = mqtot
109          etot0 = etot 
110          angtot0  = angtot
111          stot0 = stot 
112          AAM_mass_old=AAM_mass
113          AAM_vel_old=AAM_vel
114          AAM_velp_old=AAM_velp
115          AAM_velm_old=AAM_velm
116          AAM_mass_source=0.
117          AAM_vel_source=0.
118          AAM_vel_plus_source=0.
119          AAM_vel_minus_source=0.
120       END IF
121       rmsvtot=SQRT(rmsvtot/mtot) 
122       ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1.
123       mqtot=(mqtot+1E-100)/(mqtot0+1e-100)-1.
124       etot=etot/etot0-1. ; angtot=angtot/angtot0-1. ; stot=stot/stot0-1. 
125       rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo) 
126       
127       OPEN(134,file="checkconsicosa.txt",position='append')
128       WRITE(134,4000)mtot,mqtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
129       WRITE(134,*)mtot,mqtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot 
130       WRITE(134,*)"=================================================="
131       WRITE(*,4000)mtot,mqtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
132       
1334000   FORMAT(10x,'masse',5x,'advec mass',5x,'rmsdpdt',6x,'energie',3x,'enstrophie'  &
134            ,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB  '          &
135            ,2e10.3,e13.5,5e13.3/)     
136       CLOSE(134)
137
138       IF(check_type == check_detailed) THEN
139          !ang_mass = AAM_mass-AAM_mass_old
140          !WRITE(*,'(A,3E12.4)') 'AAM_mass sanity check', SUM(AAM_mass_source), ang_mass, SUM(AAM_mass_source)-ang_mass
141          !ang_vel = AAM_vel-AAM_vel_old
142          !WRITE(*,'(A,3E12.4)') 'AAM_vel  sanity check', SUM(AAM_vel_source), ang_vel, SUM(AAM_vel_source)-ang_vel
143          WRITE(*,'(A,6E12.4)') 'AAM_mass : time,old,new,dissip,dyn,phys', dt*it, AAM_mass_old, AAM_mass, &
144               AAM_mass_source(AAM_dissip), AAM_mass_source(AAM_dyn), AAM_mass_source(AAM_phys)
145          WRITE(*,'(A,6E12.4)') 'AAM_vel  : time,old,new,dissip,dyn,phys', dt*it, AAM_vel_old, AAM_vel, &
146               AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys)
147          WRITE(*,'(A,6E12.4)') 'AAM_vel+ : time,old,new,dissip,dyn,phys', dt*it, AAM_velp_old, AAM_velp, &
148               AAM_vel_plus_source(AAM_dissip), AAM_vel_plus_source(AAM_dyn), AAM_vel_plus_source(AAM_phys)
149          WRITE(*,'(A,6E12.4)') 'AAM_vel- : time,old,new,dissip,dyn,phys', dt*it, AAM_velm_old, AAM_velm, &
150               AAM_vel_minus_source(AAM_dissip), AAM_vel_minus_source(AAM_dyn), AAM_vel_minus_source(AAM_phys)
151          WRITE(*,'(A,6E12.4)') 'AAM_dyn budget mass + vel:', AAM_mass_source(AAM_dyn) + AAM_vel_source(AAM_dyn)
152          AAM_mass_old=AAM_mass
153          AAM_vel_old=AAM_vel
154          AAM_velp_old=AAM_velp
155          AAM_velm_old=AAM_velm
156          AAM_mass_source=0.
157          AAM_vel_source=0.
158          AAM_vel_minus_source=0.
159          AAM_vel_plus_source=0.
160       END IF
161    END IF
162
163  END SUBROUTINE check_conserve
164
165!------------------------ Detailed check (Lauritzen et al., 2014) -----------------------------
166
167  SUBROUTINE check_conserve_detailed(it,tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis) 
168    USE pression_mod
169    USE vorticity_mod
170    USE caldyn_gcm_mod
171    USE exner_mod 
172!    USE mpipara, ONLY : is_mpi_root, comm_icosa
173    USE omp_para, ONLY : is_master
174    INTEGER, INTENT(IN) :: tag
175    TYPE(t_field),POINTER :: f_ps(:)
176    TYPE(t_field),POINTER :: f_dps(:)
177    TYPE(t_field),POINTER :: f_ue(:)
178    TYPE(t_field),POINTER :: f_theta_rhodz(:)
179    TYPE(t_field),POINTER :: f_phis(:)
180    INTEGER, INTENT(IN) ::it
181
182    REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
183    INTEGER::ind
184    REAL(rstd) ::  etot,stot,rmsv, ang_mass, ang_vel, ang_velp, ang_velm
185   
186    IF(check_type == check_detailed) THEN
187       CALL abort_acc("!check_detailed")
188       CALL transfert_request(f_ue,req_e1_vect)
189       CALL pression(f_ps,f_p)
190       DO ind=1,ndomain 
191          IF (.NOT. assigned_domain(ind)) CYCLE
192          CALL swap_dimensions(ind)
193          CALL swap_geometry(ind)
194          p=f_p(ind) 
195          rhodz=f_rhodz(ind)
196          CALL compute_rhodz(p,rhodz) 
197       END DO
198       CALL exner(f_ps,f_p,f_pks,f_pk)
199       CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsv) 
200       
201       IF (is_master) THEN
202          ! XXX=mass,vel
203          ! AAM_XXX = value of total AAM_XXX at last check
204          ! ang_XXX = value of total AAM_XXX now
205          ! AAM_XXX_source(tag) = source of AAA_XXX due to tag=dissip,dyn,phys
206          !WRITE(*,'(A,3E12.2)') 'AAM_mass', tag, ang_mass, AAM_mass, ang_mass-AAM_mass ! FIXME
207          !WRITE(*,'(A,3E12.2)') 'AAM_vel ', tag, ang_vel, AAM_vel, ang_vel-AAM_vel ! FIXME
208          AAM_mass_source(tag) = AAM_mass_source(tag) + ang_mass - AAM_mass
209          AAM_vel_source(tag) = AAM_vel_source(tag) + ang_vel - AAM_vel
210          AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_velp - AAM_velp
211          AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_velm - AAM_velm
212          !if ((ang_vel - AAM_vel) .gt. 0.) then
213          !  AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_vel - AAM_vel
214          !else
215          !  AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_vel - AAM_vel
216          !endif
217          AAM_mass = ang_mass
218          AAM_vel = ang_vel
219          AAM_velp = ang_velp
220          AAM_velm = ang_velm
221       END IF
222    END IF
223  END SUBROUTINE check_conserve_detailed
224 
225!---------------------------------------------------------------------
226
227  SUBROUTINE global_sum(loc_sum, glob_sum)
228    USE mpi_mod
229    USE mpipara
230    USE transfert_omp_mod
231    REAL(rstd), INTENT(IN) :: loc_sum
232    REAL(rstd), INTENT(OUT) :: glob_sum
233    REAL(rstd) :: sum
234
235    CALL reduce_sum_omp(loc_sum, sum)
236    IF (using_mpi) THEN
237!$OMP MASTER   
238      CALL MPI_REDUCE(sum, glob_sum, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
239!$OMP END MASTER
240    ELSE
241      glob_sum=loc_sum
242    ENDIF
243  END SUBROUTINE global_sum
244
245  SUBROUTINE check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
246  USE mpi_mod
247  USE mpipara
248  USE transfert_omp_mod
249  USE omp_para
250  IMPLICIT NONE
251    TYPE(t_field),POINTER :: f_ps(:)
252    TYPE(t_field),POINTER :: f_dps(:)
253    REAL(rstd),POINTER :: ps(:),dps(:) 
254    REAL(rstd), INTENT(OUT) :: mtot, rmsdpdt
255
256    INTEGER :: ind,i,j,ij 
257    REAL :: mloc, rmsloc
258
259    mloc=0.0; rmsloc=0.0
260    DO ind=1,ndomain
261      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
262      CALL swap_dimensions(ind)
263      CALL swap_geometry(ind)
264      dps=f_dps(ind) 
265      ps=f_ps(ind)
266      DO j=jj_begin,jj_end
267        DO i=ii_begin,ii_end
268          ij=(j-1)*iim+i
269          IF (domain(ind)%own(i,j)) THEN
270             mloc=mloc+ps(ij)*Ai(ij)
271             rmsloc=rmsloc+dps(ij)*dps(ij)
272          ENDIF
273        ENDDO
274      ENDDO
275    ENDDO
276   
277    CALL global_sum(mloc, mtot)
278    CALL global_sum(rmsloc, rmsdpdt)
279  END SUBROUTINE check_mass_conserve
280
281  SUBROUTINE check_qmass_conserve(f_q,f_rhodz,mqtot)
282  USE mpi_mod
283  USE mpipara
284  USE transfert_omp_mod
285  USE omp_para
286  IMPLICIT NONE
287    TYPE(t_field),POINTER :: f_q(:)
288    TYPE(t_field),POINTER :: f_rhodz(:)
289    REAL(rstd),POINTER :: q(:,:,:) 
290    REAL(rstd),POINTER :: rhodz(:,:) 
291    REAL(rstd), INTENT(OUT) :: mqtot
292
293    INTEGER :: ind,i,j,ij,l 
294    REAL :: mloc
295
296    mloc=0.0
297    DO ind=1,ndomain
298      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
299      CALL swap_dimensions(ind)
300      CALL swap_geometry(ind)
301      q=f_q(ind)
302      rhodz=f_rhodz(ind)
303      DO j=jj_begin,jj_end
304        DO i=ii_begin,ii_end
305          ij=(j-1)*iim+i
306          IF (domain(ind)%own(i,j)) THEN
307            DO l=1, llm
308             mloc=mloc + sum( q(ij,l,:)*Ai(ij)*rhodz(ij,l) )
309            END DO
310          ENDIF
311        ENDDO
312      ENDDO
313    ENDDO
314   
315    CALL global_sum(mloc, mqtot)
316  END SUBROUTINE check_qmass_conserve
317
318!---------------------------------------------------------------------
319
320  SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
321                      stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot) 
322  USE pression_mod
323  USE vorticity_mod
324    TYPE(t_field), POINTER :: f_ue(:)
325    TYPE(t_field), POINTER :: f_theta_rhodz(:)
326    TYPE(t_field), POINTER :: f_phis(:)
327    REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot
328    REAL(rstd), POINTER :: ue(:,:)
329    REAL(rstd), POINTER :: p(:,:) 
330    REAL(rstd), POINTER :: theta_rhodz(:,:,:) 
331    REAL(rstd), POINTER :: pk(:,:) 
332    REAL(rstd), POINTER :: phis(:) 
333    REAL(rstd), POINTER :: rhodz(:,:) 
334    REAL(rstd) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv
335    INTEGER :: ind
336
337    e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0
338    AAM_velp = 0 ; AAM_velm = 0
339
340    DO ind=1,ndomain 
341      IF (.NOT. assigned_domain(ind)) CYCLE
342      CALL swap_dimensions(ind)
343      CALL swap_geometry(ind)
344      ue=f_ue(ind) 
345      p=f_p(ind)
346      rhodz=f_rhodz(ind) 
347      theta_rhodz=f_theta_rhodz(ind) 
348      pk=f_pk(ind) 
349      phis=f_phis(ind) 
350      CALL compute_energy(ind,ue,p,rhodz,theta_rhodz(:,:,1),pk,phis, &
351           e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv)
352    END DO
353
354    CALL global_sum(e, etot)
355    CALL global_sum(s, stot)
356    CALL global_sum(rmsv, rmsvtot)
357    CALL global_sum(AAM_mass, AAM_mass_tot)
358    CALL global_sum(AAM_vel, AAM_vel_tot)
359    CALL global_sum(AAM_velp, AAM_velp_tot)
360    CALL global_sum(AAM_velm, AAM_velm_tot)
361  END SUBROUTINE check_energy
362
363  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv)
364    USE disvert_mod
365    USE wind_mod
366    USE omp_para
367    INTEGER,INTENT(IN)::ind 
368    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
369    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
370    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
371    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
372    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
373    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
374    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv
375   
376    REAL(rstd) :: ucenter(iim*jjm,llm,3)
377    REAL(rstd) :: ulon(iim*jjm,llm)
378    REAL(rstd) :: ulat(iim*jjm,llm)
379   
380    REAL(rstd) :: mass_ij, theta_ij, KE_ij, rad,radomeg,lat
381    INTEGER :: i,j,ij,l
382   
383    CALL compute_wind_centered(u,ucenter)
384    CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
385    rad = radius 
386    radomeg = rad*omega
387   
388    DO l = ll_begin,ll_end
389       DO j=jj_begin,jj_end
390          DO i=ii_begin,ii_end
391             ij=(j-1)*iim+i
392             IF (domain(ind)%own(i,j)) THEN
393                lat = lat_i(ij) 
394                mass_ij = rhodz(ij,l)*Ai(ij) 
395                theta_ij = theta_rhodz(ij,l)/rhodz(ij,l)
396                KE_ij = 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
397                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
398                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
399                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
400                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
401                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
402                rmsv = rmsv + 2*mass_ij*KE_ij
403                s    = s + Ai(ij)*theta_rhodz(ij,l)
404                e    = e + mass_ij*(phis(ij)+theta_ij*pk(ij,l)+KE_ij)
405                AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat)
406                AAM_vel  = AAM_vel  + rad*cos(lat)*mass_ij*ulon(ij,l)
407                if (ulon(ij,l) > 0.) then
408                  AAM_velp  = AAM_velp  + rad*cos(lat)*mass_ij*ulon(ij,l)
409                else
410                  AAM_velm  = AAM_velm  + rad*cos(lat)*mass_ij*ulon(ij,l)
411                endif
412             END IF
413          END DO
414       END DO
415    END DO
416   
417  END SUBROUTINE compute_energy
418
419!---------------------------------------------------------------------
420 
421  SUBROUTINE check_PV(ztot)
422    REAL(rstd),INTENT(OUT) :: ztot
423    REAL(rstd), POINTER :: vort(:,:)
424    REAL(rstd), POINTER :: rhodz(:,:) 
425    INTEGER :: ind
426    REAL(rstd) :: z
427   
428    z=0
429    DO ind=1,ndomain 
430       IF (.NOT. assigned_domain(ind)) CYCLE
431       CALL swap_dimensions(ind)
432       CALL swap_geometry(ind)
433       vort=f_vort(ind)
434       rhodz=f_rhodz(ind) 
435       CALL compute_PV(vort,rhodz,z) 
436    ENDDO
437    CALL global_sum(z, ztot)
438  END SUBROUTINE check_PV 
439
440  SUBROUTINE compute_PV(vort,rhodz,z) 
441  USE disvert_mod
442  USE omp_para
443    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
444    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
445    REAL(rstd),INTENT(INOUT) :: z
446    REAL(rstd)::qv1,qv2 
447    REAL(rstd)::hv1,hv2 
448    INTEGER :: i,j,ij,l
449
450    hv1 = 0.0 ; hv2 = 0.0 
451   
452    DO l = ll_begin,ll_end
453      DO j=jj_begin+1,jj_end-1
454        DO i=ii_begin+1,ii_end-1
455          ij=(j-1)*iim+i
456          hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
457            + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
458            + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
459   
460          qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
461         
462          hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
463            + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
464            + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
465                     
466          qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
467
468          z = z + (qv1*qv1*hv1 + qv2*qv2*hv2) 
469
470        ENDDO
471      ENDDO
472     
473    ENDDO
474
475  END SUBROUTINE compute_PV
476!---------------------------------------------------------------------
477
478   SUBROUTINE compute_rhodz(p,rhodz) 
479     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
480     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
481     REAL(rstd) ::mass(iim*jjm,llm)
482     INTEGER :: i,j,ij,l 
483       
484     DO l = 1, llm
485      DO j=jj_begin,jj_end
486       DO i=ii_begin,ii_end
487         ij=(j-1)*iim+i
488         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
489         rhodz(ij,l) = mass(ij,l) / Ai(ij)
490       ENDDO
491      ENDDO
492     ENDDO
493
494   END SUBROUTINE compute_rhodz
495!====================================================================
496
497
498 END MODULE check_conserve_mod
Note: See TracBrowser for help on using the repository browser.