Changeset 326 for codes/icosagcm/trunk
- Timestamp:
- 02/02/15 01:12:53 (9 years ago)
- Location:
- codes/icosagcm/trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/trunk/src/check_conserve.f90
r295 r326 7 7 TYPE(t_field),POINTER,SAVE :: f_vort(:) 8 8 TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 9 10 PUBLIC init_check_conserve, check_conserve 9 10 INTEGER, PARAMETER :: check_basic=1, check_detailed=2 11 INTEGER, SAVE :: check_type 12 PUBLIC init_check_conserve, check_conserve, check_conserve_detailed 11 13 REAL(rstd),SAVE:: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0 12 !$OMP THREADPRIVATE( mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)14 !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0) 13 15 14 16 CONTAINS … … 17 19 18 20 SUBROUTINE init_check_conserve 19 USE icosa 20 IMPLICIT NONE 21 USE icosa 22 USE getin_mod 23 USE omp_para, ONLY : is_master 24 IMPLICIT NONE 25 CHARACTER(LEN=255) :: check_type_str 21 26 CALL allocate_field(f_pk,field_t,type_real,llm) 22 27 CALL allocate_field(f_p,field_t,type_real,llm+1) … … 24 29 CALL allocate_field(f_rhodz,field_t,type_real,llm) 25 30 CALL allocate_field(f_vort,field_z,type_real,llm) 31 32 check_type_str='basic' 33 CALL getin("check_conservation",check_type_str) 34 SELECT CASE(TRIM(check_type_str)) 35 CASE('basic') 36 check_type=check_basic 37 CASE('detailed') 38 check_type=check_detailed 39 CASE DEFAULT 40 IF(is_master) PRINT*,'Bad selector <', TRIM(check_type_str), & 41 '> for variable check_conservation : options are <basic>, <detailed>' 42 STOP 43 END SELECT 26 44 END SUBROUTINE init_check_conserve 45 46 !--------------------------------- Basic check -------------------------------- 27 47 28 48 SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) … … 40 60 TYPE(t_field),POINTER :: f_theta_rhodz(:) 41 61 TYPE(t_field),POINTER :: f_phis(:) 42 INTEGER ::it62 INTEGER, INTENT(IN) :: it 43 63 44 64 REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 45 65 INTEGER::ind,ierr 46 REAL(rstd) :: mtot, rmsdpdt47 REAL(rstd) :: etot, stot, ang tot, rmsvtot, ztot66 REAL(rstd) :: mtot, angtot, rmsdpdt 67 REAL(rstd) :: etot, stot, ang_mass, ang_vit, rmsvtot, ztot 48 68 49 69 CALL transfert_request(f_ue,req_e1_vect) … … 63 83 CALL check_PV(ztot) 64 84 CALL exner(f_ps,f_p,f_pks,f_pk) 65 CALL check_EN(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot) 85 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vit, rmsvtot) 86 angtot = ang_mass + ang_vit 66 87 67 88 IF (is_master) THEN … … 70 91 mtot0 = mtot 71 92 etot0 = etot 72 angtot0 = angtot 93 angtot0 = angtot 73 94 stot0 = stot 74 95 END IF … … 92 113 END IF 93 114 END SUBROUTINE check_conserve 94 115 116 !------------------------ Detailed check (Lauritzen et al., 2014) ----------------------------- 117 118 SUBROUTINE check_conserve_detailed(tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 119 USE icosa 120 USE pression_mod 121 USE vorticity_mod 122 USE caldyn_gcm_mod 123 USE exner_mod 124 USE mpipara, ONLY : is_mpi_root, comm_icosa 125 IMPLICIT NONE 126 CHARACTER(LEN=*), INTENT(IN) :: tag 127 TYPE(t_field),POINTER :: f_ps(:) 128 TYPE(t_field),POINTER :: f_dps(:) 129 TYPE(t_field),POINTER :: f_ue(:) 130 TYPE(t_field),POINTER :: f_theta_rhodz(:) 131 TYPE(t_field),POINTER :: f_phis(:) 132 INTEGER, INTENT(IN) ::it 133 134 REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 135 INTEGER::ind,ierr 136 REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, AAM_mass, AAM_vel 137 138 IF(check_type == check_detailed) THEN 139 etot=0.0; stot=0.0; rmsv=0.0 140 AAM_vel=0. ; AAM_vel=0. 141 ztot = 0.0 142 143 CALL transfert_request(f_ue,req_e1_vect) 144 CALL pression(f_ps,f_p) 145 DO ind=1,ndomain 146 IF (.NOT. assigned_domain(ind)) CYCLE 147 CALL swap_dimensions(ind) 148 CALL swap_geometry(ind) 149 p=f_p(ind) 150 rhodz=f_rhodz(ind) 151 CALL compute_rhodz(p,rhodz) 152 END DO 153 CALL exner(f_ps,f_p,f_pks,f_pk) 154 CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, AAM_mass, AAM_vel, rmsv) 155 156 IF (is_mpi_root) THEN 157 !$OMP MASTER 158 IF ( it == itau0 ) Then 159 ztot0 = ztot 160 mtot0 = mtot 161 etot0 = etot 162 angtot0 = AAM_mass + AAM_vel 163 stot0 = stot 164 END IF 165 166 PRINT *, TRIM(tag), AAM_mass, AAM_vel 167 !$OMP END MASTER 168 END IF 169 END IF 170 END SUBROUTINE check_conserve_detailed 171 95 172 !--------------------------------------------------------------------- 173 174 SUBROUTINE global_sum(loc_sum, glob_sum) 175 USE mpi_mod 176 USE mpipara 177 USE transfert_omp_mod 178 IMPLICIT NONE 179 REAL(rstd), INTENT(IN) :: loc_sum 180 REAL(rstd), INTENT(OUT) :: glob_sum 181 REAL(rstd) :: sum 182 183 CALL reduce_sum_omp(loc_sum, sum) 184 IF (using_mpi) THEN 185 !$OMP MASTER 186 CALL MPI_REDUCE(sum, glob_sum, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 187 !$OMP END MASTER 188 ELSE 189 glob_sum=loc_sum 190 ENDIF 191 END SUBROUTINE global_sum 96 192 97 193 SUBROUTINE check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt) … … 128 224 ENDDO 129 225 ENDDO 130 131 IF (using_mpi) THEN 132 CALL reduce_sum_omp(mloc, mloc_mpi) 133 CALL reduce_sum_omp(rmsloc, rmsloc_mpi) 134 !$OMP MASTER 135 CALL MPI_REDUCE(mloc_mpi, mtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 136 CALL MPI_REDUCE(rmsloc_mpi, rmsdpdt, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 137 !$OMP END MASTER 138 ELSE 139 mtot=mloc 140 rmsdpdt=rmsloc 141 ENDIF 142 226 227 CALL global_sum(mloc, mtot) 228 CALL global_sum(rmsloc, rmsdpdt) 143 229 END SUBROUTINE check_mass_conserve 144 230 145 231 !--------------------------------------------------------------------- 146 232 147 SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot) 233 SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, & 234 stot, AAM_mass_tot, AAM_vel_tot, rmsvtot) 148 235 USE icosa 149 236 USE pression_mod 150 237 USE vorticity_mod 151 USE mpi_mod152 USE mpipara153 USE transfert_omp_mod154 238 IMPLICIT NONE 155 239 TYPE(t_field), POINTER :: f_ue(:) 156 240 TYPE(t_field), POINTER :: f_theta_rhodz(:) 157 241 TYPE(t_field), POINTER :: f_phis(:) 158 REAL(rstd),INTENT(OUT) :: etot, stot, angtot, rmsvtot242 REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot 159 243 REAL(rstd), POINTER :: ue(:,:) 160 244 REAL(rstd), POINTER :: p(:,:) … … 163 247 REAL(rstd), POINTER :: phis(:) 164 248 REAL(rstd), POINTER :: rhodz(:,:) 165 REAL(rstd) :: e, s, ang, rmsv 166 REAL(rstd) :: e_mpi, s_mpi, ang_mpi, rmsv_mpi 249 REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv 167 250 INTEGER :: ind 168 251 169 e = 0 ; s = 0 ; ang =0 ; rmsv = 0252 e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0 170 253 171 254 DO ind=1,ndomain … … 179 262 pk=f_pk(ind) 180 263 phis=f_phis(ind) 181 CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv) 264 CALL compute_energy(ind,ue,p,rhodz,theta_rhodz,pk,phis, & 265 e, s, AAM_mass, AAM_vel, rmsv) 182 266 END DO 183 267 184 IF (using_mpi) THEN 185 CALL reduce_sum_omp(e, e_mpi) 186 CALL reduce_sum_omp(s, s_mpi) 187 CALL reduce_sum_omp(ang, ang_mpi) 188 CALL reduce_sum_omp(rmsv, rmsv_mpi) 189 !$OMP MASTER 190 CALL MPI_REDUCE(e_mpi, etot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 191 CALL MPI_REDUCE(s_mpi, stot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 192 CALL MPI_REDUCE(ang_mpi, angtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 193 CALL MPI_REDUCE(rmsv_mpi, rmsvtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr) 194 !$OMP END MASTER 195 ELSE 196 etot=e 197 stot=s 198 angtot=ang 199 rmsvtot=rmsv 200 ENDIF 201 202 END SUBROUTINE check_en 203 204 SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv) 205 USE icosa 206 USE disvert_mod 207 USE wind_mod 208 USE omp_para 209 IMPLICIT NONE 268 CALL global_sum(e, etot) 269 CALL global_sum(s, stot) 270 CALL global_sum(rmsv, rmsvtot) 271 CALL global_sum(AAM_mass, AAM_mass_tot) 272 CALL global_sum(AAM_vel, AAM_vel_tot) 273 END SUBROUTINE check_energy 274 275 SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv) 276 USE icosa 277 USE disvert_mod 278 USE wind_mod 279 USE omp_para 280 IMPLICIT NONE 210 281 INTEGER,INTENT(IN)::ind 211 282 REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) … … 215 286 REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm) 216 287 REAL(rstd),INTENT(IN) :: phis(iim*jjm) 217 REAL(rstd),INTENT(INOUT) :: e 218 REAL(rstd),INTENT(INOUT) :: s 219 REAL(rstd),INTENT(INOUT) :: ang 220 REAL(rstd),INTENT(INOUT) :: rmsv 221 222 REAL(rstd):: theta(iim*jjm,llm) 223 REAL(rstd)::KE(iim*jjm,llm) 288 REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv 289 224 290 REAL(rstd) :: ucenter(iim*jjm,3,llm) 225 291 REAL(rstd) :: ulon(iim*jjm,llm) 226 292 REAL(rstd) :: ulat(iim*jjm,llm) 227 228 REAL(rstd)::masse(iim*jjm,llm) 229 REAL(rstd)::rad,radomeg,lat,lon 230 REAL(rstd) ::xyz(3) 231 INTEGER :: i,j,ij,l,ij2 232 233 234 DO l = ll_begin, ll_end 235 DO j=jj_begin,jj_end 236 DO i=ii_begin,ii_end 237 ij=(j-1)*iim+i 238 masse(ij,l) = rhodz(ij,l)*Ai(ij) 239 theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 240 IF (domain(ind)%own(i,j)) THEN 241 s = s + Ai(ij)*theta_rhodz(ij,l) 242 END IF 243 END DO 244 END DO 245 END DO 246 247 DO l = ll_begin,ll_end 248 DO j=jj_begin,jj_end 249 DO i=ii_begin,ii_end 250 ij=(j-1)*iim+i 251 IF (domain(ind)%own(i,j)) THEN 252 KE(ij,l)= 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & 253 le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 + & 254 le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 + & 255 le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 + & 256 le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 257 le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 258 rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 259 END IF 260 END DO 261 END DO 262 END DO 263 264 DO l = ll_begin,ll_end 265 DO j=jj_begin,jj_end 266 DO i=ii_begin,ii_end 267 ij=(j-1)*iim+i 268 IF (domain(ind)%own(i,j)) THEN 269 e = e + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l)) 270 END IF 271 END DO 272 END DO 273 END DO 274 275 !------------------------------ ANGULAR VELOCITY 293 294 REAL(rstd) :: mass_ij, theta_ij, KE_ij, rad,radomeg,lat 295 INTEGER :: i,j,ij,l 296 276 297 CALL compute_wind_centered(u,ucenter) 277 298 CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat) 278 299 rad = radius 279 300 radomeg = rad*omega 280 301 281 302 DO l = ll_begin,ll_end 282 DO j=jj_begin,jj_end 283 DO i=ii_begin,ii_end 284 ij=(j-1)*iim+i 285 IF (domain(ind)%own(i,j)) THEN 286 lat=lat_i(ij) 287 ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat)) 288 END IF 289 END DO 290 END DO 291 END DO 292 293 END SUBROUTINE compute_en 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 lat = lat_i(ij) 308 mass_ij = rhodz(ij,l)*Ai(ij) 309 theta_ij = theta_rhodz(ij,l)/rhodz(ij,l) 310 KE_ij = 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 + & 311 le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 + & 312 le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 + & 313 le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 + & 314 le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 + & 315 le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 316 rmsv = rmsv + 2*mass_ij*KE_ij 317 s = s + Ai(ij)*theta_rhodz(ij,l) 318 e = e + mass_ij*(phis(ij)+theta_ij*pk(ij,l)+KE_ij) 319 AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat) 320 AAM_vel = AAM_vel + rad*cos(lat)*mass_ij*ulon(ij,l) 321 END IF 322 END DO 323 END DO 324 END DO 325 326 END SUBROUTINE compute_energy 294 327 295 328 !--------------------------------------------------------------------- -
codes/icosagcm/trunk/src/timeloop_gcm.f90
r295 r326 211 211 212 212 DO it=itau0+1,itau0+itaumax 213 213 214 CALL check_conserve_detailed('detailed_budget 0', & 215 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 216 214 217 IF (xios_output) CALL xios_update_calendar(it) 215 218 IF (it==itau0+1 .OR. MOD(it,itau_sync)==0) THEN … … 256 259 END SELECT 257 260 END DO 261 262 CALL check_conserve_detailed('detailed_budget 1', & 263 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 258 264 259 265 IF (MOD(it,itau_dissip)==0) THEN … … 281 287 END IF 282 288 289 CALL check_conserve_detailed('detailed_budget 2', & 290 f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 291 283 292 IF(MOD(it,itau_adv)==0) THEN 284 293
Note: See TracChangeset
for help on using the changeset viewer.