Changeset 649
- Timestamp:
- 12/21/17 22:55:38 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/diagnostics/check_conserve.f90
r607 r649 14 14 15 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, & 16 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) 17 19 REAL(rstd),SAVE :: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0 18 20 !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0) … … 66 68 INTEGER :: ind,ierr 67 69 REAL(rstd) :: mtot, angtot, rmsdpdt 68 REAL(rstd) :: etot, stot, ang_mass, ang_vel, rmsvtot, ztot70 REAL(rstd) :: etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot, ztot 69 71 70 72 CALL transfert_request(f_ue,req_e1_vect) … … 84 86 CALL check_PV(ztot) 85 87 CALL exner(f_ps,f_p,f_pks,f_pk) 86 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsvtot)88 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot) 87 89 88 90 IF (is_master) THEN ! is_master = is_omp_master && is_mpi_master … … 90 92 AAM_mass = ang_mass 91 93 AAM_vel = ang_vel 94 AAM_velp = ang_velp 95 AAM_velm = ang_velm 92 96 angtot = ang_mass + ang_vel 93 97 IF ( it == itau0 ) THEN … … 99 103 AAM_mass_old=AAM_mass 100 104 AAM_vel_old=AAM_vel 105 AAM_velp_old=AAM_velp 106 AAM_velm_old=AAM_velm 101 107 AAM_mass_source=0. 102 108 AAM_vel_source=0. 109 AAM_vel_plus_source=0. 110 AAM_vel_minus_source=0. 103 111 END IF 104 112 rmsvtot=SQRT(rmsvtot/mtot) … … 127 135 WRITE(*,'(A,6E12.4)') 'AAM_vel : time,old,new,dissip,dyn,phys', dt*it, AAM_vel_old, AAM_vel, & 128 136 AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys) 129 WRITE(*,'(A,6E12.4)') 'AAM_ tot : time,old,new,dissip,dyn,phys', dt*it, &130 AAM_vel_ old+AAM_mass_old, AAM_vel+AAM_mass, &131 AAM_vel_source(AAM_dissip)+AAM_mass_source(AAM_dissip), &132 AAM_vel_ source(AAM_dyn) +AAM_mass_source(AAM_dyn), &133 AAM_vel_source(AAM_phys) +AAM_mass_source(AAM_phys)137 WRITE(*,'(A,6E12.4)') 'AAM_vel+ : time,old,new,dissip,dyn,phys', dt*it, AAM_velp_old, AAM_velp, & 138 AAM_vel_plus_source(AAM_dissip), AAM_vel_plus_source(AAM_dyn), AAM_vel_plus_source(AAM_phys) 139 WRITE(*,'(A,6E12.4)') 'AAM_vel- : time,old,new,dissip,dyn,phys', dt*it, AAM_velm_old, AAM_velm, & 140 AAM_vel_minus_source(AAM_dissip), AAM_vel_minus_source(AAM_dyn), AAM_vel_minus_source(AAM_phys) 141 WRITE(*,'(A,6E12.4)') 'AAM_dyn budget mass + vel:', AAM_mass_source(AAM_dyn) + AAM_vel_source(AAM_dyn) 134 142 AAM_mass_old=AAM_mass 135 143 AAM_vel_old=AAM_vel 144 AAM_velp_old=AAM_velp 145 AAM_velm_old=AAM_velm 136 146 AAM_mass_source=0. 137 147 AAM_vel_source=0. 148 AAM_vel_minus_source=0. 149 AAM_vel_plus_source=0. 138 150 END IF 139 151 END IF … … 160 172 REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 161 173 INTEGER::ind,ierr 162 REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel 174 REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel, ang_velp, ang_velm 163 175 164 176 IF(check_type == check_detailed) THEN … … 175 187 END DO 176 188 CALL exner(f_ps,f_p,f_pks,f_pk) 177 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsv)189 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsv) 178 190 179 191 IF (is_master) THEN … … 186 198 AAM_mass_source(tag) = AAM_mass_source(tag) + ang_mass - AAM_mass 187 199 AAM_vel_source(tag) = AAM_vel_source(tag) + ang_vel - AAM_vel 200 AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_velp - AAM_velp 201 AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_velm - AAM_velm 202 !if ((ang_vel - AAM_vel) .gt. 0.) then 203 ! AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_vel - AAM_vel 204 !else 205 ! AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_vel - AAM_vel 206 !endif 188 207 AAM_mass = ang_mass 189 208 AAM_vel = ang_vel 209 AAM_velp = ang_velp 210 AAM_velm = ang_velm 190 211 END IF 191 212 END IF … … 252 273 253 274 SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, & 254 stot, AAM_mass_tot, AAM_vel_tot, rmsvtot)275 stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot) 255 276 USE pression_mod 256 277 USE vorticity_mod … … 258 279 TYPE(t_field), POINTER :: f_theta_rhodz(:) 259 280 TYPE(t_field), POINTER :: f_phis(:) 260 REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot281 REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot 261 282 REAL(rstd), POINTER :: ue(:,:) 262 283 REAL(rstd), POINTER :: p(:,:) … … 265 286 REAL(rstd), POINTER :: phis(:) 266 287 REAL(rstd), POINTER :: rhodz(:,:) 267 REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv288 REAL(rstd) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv 268 289 INTEGER :: ind 269 290 270 291 e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0 292 AAM_velp = 0 ; AAM_velm = 0 271 293 272 294 DO ind=1,ndomain … … 281 303 phis=f_phis(ind) 282 304 CALL compute_energy(ind,ue,p,rhodz,theta_rhodz(:,:,1),pk,phis, & 283 e, s, AAM_mass, AAM_vel, rmsv)305 e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv) 284 306 END DO 285 307 … … 289 311 CALL global_sum(AAM_mass, AAM_mass_tot) 290 312 CALL global_sum(AAM_vel, AAM_vel_tot) 313 CALL global_sum(AAM_velp, AAM_velp_tot) 314 CALL global_sum(AAM_velm, AAM_velm_tot) 291 315 END SUBROUTINE check_energy 292 316 293 SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv)317 SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv) 294 318 USE disvert_mod 295 319 USE wind_mod … … 302 326 REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm) 303 327 REAL(rstd),INTENT(IN) :: phis(iim*jjm) 304 REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv328 REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv 305 329 306 330 REAL(rstd) :: ucenter(iim*jjm,llm,3) … … 335 359 AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat) 336 360 AAM_vel = AAM_vel + rad*cos(lat)*mass_ij*ulon(ij,l) 361 if (ulon(ij,l) > 0.) then 362 AAM_velp = AAM_velp + rad*cos(lat)*mass_ij*ulon(ij,l) 363 else 364 AAM_velm = AAM_velm + rad*cos(lat)*mass_ij*ulon(ij,l) 365 endif 337 366 END IF 338 367 END DO
Note: See TracChangeset
for help on using the changeset viewer.