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

Last change on this file since 930 was 902, checked in by adurocher, 5 years ago

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