Changeset 1569 for trunk/LMDZ.MARS/libf/aeronomars/photochemistry_asis.F90
- Timestamp:
- Jun 29, 2016, 6:08:29 PM (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/aeronomars/photochemistry_asis.F90
r1503 r1569 85 85 integer,parameter :: i_n2 = 17 86 86 87 real :: stimestep ! standard timestep for the chemistry (s) 88 real :: ctimestep ! real timestep for the chemistry (s) 87 real :: ctimestep ! standard timestep for the chemistry (s) 89 88 real :: dt_guess ! first-guess timestep (s) 90 89 real :: dt_corrected ! corrected timestep (s) 91 real :: dt_min = 1. ! minimum allowed timestep (s)92 real :: dtg ! correction factor for the timestep (s)93 90 real :: j(nlayer,nd) ! interpolated photolysis rates (s-1) 94 91 real :: time ! internal time (between 0 and ptimestep, in s) … … 115 112 116 113 real (kind = 8), dimension(nesp) :: prod, loss 117 118 ! curvatures119 120 real :: ratio, curv, e, e1, e2, e3121 114 122 115 !=================================================================== … … 169 162 170 163 !=================================================================== 171 ! stimestep : standard chemical timestep (s) 172 ! ctimestep : real chemical timestep (s), 173 ! taking into account the physical timestep 174 !=================================================================== 175 176 stimestep = 600. ! standard value : 10 mn 177 178 phychemrat = nint(ptimestep/stimestep) 164 ! ctimestep : standard chemical timestep (s), defined as 165 ! the fraction phychemrat of the physical timestep 166 !=================================================================== 167 179 168 phychemrat = 1 180 169 181 170 ctimestep = ptimestep/real(phychemrat) 182 171 183 !print*, "stimestep = ", stimestep184 172 !print*, "ptimestep = ", ptimestep 185 173 !print*, "phychemrat = ", phychemrat … … 210 198 call fill_matrix(ilev, mat1, prod, loss, c, nesp, nlayer, v_phot, v_3, v_4) 211 199 212 ! first guess : form the matrix identity + mat*dt_guess 213 214 mat(:,:) = mat1(:,:)*dt_guess 200 ! adaptative evaluation of the sub time step 201 202 call define_dt(nesp, dt_corrected, dt_guess, ctimestep, cold(:), c(ilev,:), & 203 mat1, prod, loss) 204 205 if (time + dt_corrected > ptimestep) then 206 dt_corrected = ptimestep - time 207 end if 208 209 ! if (dt_corrected /= dt_guess) then ! the timestep has been modified 210 211 ! form the matrix identity + mat*dt_corrected 212 213 mat(:,:) = mat1(:,:)*dt_corrected 215 214 do iesp = 1,nesp 216 215 mat(iesp,iesp) = 1. + mat(iesp,iesp) 217 216 end do 218 217 219 ! first-guess: solve linear system218 ! solve the linear system M*Cn+1 = Cn (RHS in cnew, then replaced by solution) 220 219 221 220 cnew(:) = c(ilev,:) … … 227 226 stop 228 227 #endif 228 229 ! end if 230 229 231 ! eliminate small values 230 232 … … 233 235 end where 234 236 235 ! curvature of the solution and weighted error norm for O 236 ! 237 ! cold : concentrations at time t - dt 238 ! c : concentrations at time t 239 ! cnew : concentrations guessed at time t + dt 240 241 ! ratio = dtold/dtnew 242 ! here ratio = 1 since the timestep for the guess is equal to the previous timestep 243 244 ratio = 1. 245 246 curv = 2.*(ratio*cnew(i_o) - (1. + ratio)*c(ilev,i_o) + cold(i_o)) & 247 /(ratio*(1. + ratio)) 248 e1 = (curv/(cnew(i_o) + cnew(i_o3)))*100. 249 e1 = abs(e1) 250 251 ! curvature of the solution and weighted error norm for HO2 252 253 curv = 2.*(ratio*cnew(i_ho2) - (1. + ratio)*c(ilev,i_ho2) + cold(i_ho2)) & 254 /(ratio*(1. + ratio)) 255 e2 = (curv/(cnew(i_h) + cnew(i_oh) + cnew(i_ho2) + 2.*cnew(i_h2o2)))*100. 256 e2 = abs(e2) 257 258 ! curvature of the solution and weighted error norm for NO2 259 260 ! curv = 2.*(ratio*cnew(i_no2) - (1. + ratio)*c(ilev,i_no2) + cold(i_no2)) & 261 ! /(ratio*(1. + ratio)) 262 ! e3 = (curv/(cnew(i_no) + cnew(i_no2)))*100. 263 ! e3 = abs(e3) 264 e3 = 0. 265 266 ! timestep correction 267 268 e = max(e1, e2) 269 e = max(e, e3) 270 e = max(e, 0.1) 271 dtg = max(0.001, min(2.5,0.8/sqrt(e))) 272 ! dtg = 1. ! uncomment this line to turn off the variable timestep scheme 273 dt_corrected = max(dt_min,dtg*dt_guess) ! minimal timestep: dt_min 274 dt_corrected = min(dt_corrected,ctimestep) ! maximal timestep: ctimestep 275 276 if (time + dt_corrected > ptimestep) then 277 dt_corrected = ptimestep - time 278 end if 279 280 if (dt_corrected /= dt_guess) then ! the timestep has been modified 281 282 ! form the matrix identity + mat*dt_corrected 283 284 mat(:,:) = mat1(:,:)*dt_corrected 285 do iesp = 1,nesp 286 mat(iesp,iesp) = 1. + mat(iesp,iesp) 287 end do 288 289 ! solve linear system 290 291 cnew(:) = c(ilev,:) 237 ! update concentrations 238 239 cold(:) = c(ilev,:) 240 c(ilev,:) = cnew(:) 241 cnew(:) = 0. 242 243 ! increment internal time 244 245 time = time + dt_corrected 246 dt_guess = dt_corrected ! first-guess timestep for next iteration 247 248 end do ! while (time < ptimestep) 249 250 end do ! ilev 251 252 !=================================================================== 253 ! save chemical species for the gcm 254 !=================================================================== 255 256 call chimtogcm(nlayer, nq, zycol, lswitch, nesp, & 257 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 258 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & 259 i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 260 261 contains 262 263 !================================================================ 264 265 subroutine define_dt(nesp, dtnew, dtold, ctimestep, cold, ccur, mat1, prod, loss) 266 267 !================================================================ 268 ! iterative evaluation of the appropriate time step dtnew 269 ! according to curvature criterion based on 270 ! e = 2 Rtol [r Cn+1 -(1-r) Cn + Cn-1 ]/[(1+r) Cn] 271 ! with r = (tn - tn-1)/(tn+1 - tn) 272 !================================================================ 273 274 implicit none 275 276 ! input 277 278 integer :: nesp ! number of species in the chemistry 279 280 real :: dtold, ctimestep 281 real (kind = 8), dimension(nesp) :: cold, ccur 282 real (kind = 8), dimension(nesp,nesp) :: mat1 283 real (kind = 8), dimension(nesp) :: prod, loss 284 285 ! output 286 287 real :: dtnew 288 289 ! local 290 291 real (kind = 8), dimension(nesp) :: cnew 292 real (kind = 8), dimension(nesp,nesp) :: mat 293 real (kind = 8) :: ratio, e, es, coef 294 295 integer :: code, iesp, iter 296 integer, dimension(nesp) :: indx 297 298 real :: dttest 299 300 ! parameters 301 302 real (kind = 8), parameter :: dtmin = 10. ! minimum time step (s) 303 real (kind = 8), parameter :: atol = 1.e4 ! recommended value : 1.e4 304 real (kind = 8), parameter :: rtol = 1./0.1 ! 1/rtol recommended value : 0.1-0.02 305 integer, parameter :: niter = 3 ! number of iterations 306 real (kind = 8), parameter :: coefmax = 2. 307 real (kind = 8), parameter :: coefmin = 0.1 308 logical :: fast_guess = .true. 309 310 dttest = dtold ! dttest = dtold = dt_guess 311 312 do iter = 1,niter 313 314 if (fast_guess) then 315 316 ! first guess : fast implicit method 317 318 do iesp = 1, nesp 319 cnew(iesp) = (ccur(iesp) + prod(iesp)*dttest)/(1. + loss(iesp)*dttest) 320 end do 321 322 else 323 324 ! first guess : form the matrix identity + mat*dt_guess 325 326 mat(:,:) = mat1(:,:)*dttest 327 do iesp = 1,nesp 328 mat(iesp,iesp) = 1. + mat(iesp,iesp) 329 end do 330 331 ! form right-hand side (RHS) of the system 332 333 cnew(:) = ccur(:) 334 335 ! solve the linear system M*Cn+1 = Cn (RHS in cnew, then replaced by solution) 292 336 293 337 #ifdef LAPACK … … 297 341 stop 298 342 #endif 343 344 end if 345 346 ! ratio old/new subtimestep 347 348 ratio = dtold/dttest 349 350 ! e : local error indicatocitr 351 352 e = 0. 353 354 do iesp = 1,nesp 355 es = 2.*abs((ratio*cnew(iesp) - (1. + ratio)*ccur(iesp) + cold(iesp)) & 356 /(1. + ratio)/max(ccur(iesp),atol)) 357 358 if (es > e) then 359 e = es 299 360 end if 300 301 ! eliminate small values 302 303 where (cnew(:)/dens(ilev) < 1.e-30) 304 cnew(:) = 0. 305 end where 306 307 cold(:) = c(ilev,:) 308 c(ilev,:) = cnew(:) 309 cnew(:) = 0. 310 311 ! increment internal time 312 313 time = time + dt_corrected 314 dt_guess = dt_corrected ! first-guess timestep for next iteration 315 316 end do ! while (time < ptimestep) 317 318 end do ! ilev 319 320 !=================================================================== 321 ! save chemical species for the gcm 322 !=================================================================== 323 324 call chimtogcm(nlayer, nq, zycol, lswitch, nesp, & 325 i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h, & 326 i_h2, i_oh, i_ho2, i_h2o2, i_h2o, & 327 i_n, i_n2d, i_no, i_no2, i_n2, dens, c) 328 329 contains 361 end do 362 e = rtol*e 363 364 ! timestep correction 365 366 coef = max(coefmin, min(coefmax,0.8/sqrt(e))) 367 368 dttest = max(dtmin,dttest*coef) 369 dttest = min(ctimestep,dttest) 370 371 end do ! iter 372 373 ! new timestep 374 375 dtnew = dttest 376 377 end subroutine define_dt 330 378 331 379 !====================================================================== … … 973 1021 integer :: iphot,i3,i4 974 1022 975 real(kind = jprb) :: eps, eps_4 ! implicit/explicit coefficient1023 real(kind = 8) :: eps, eps_4 ! implicit/explicit coefficient 976 1024 977 1025 ! initialisations … … 1031 1079 1032 1080 eps_4 = abs(c(ilev,ind_4_2))/(abs(c(ilev,ind_4_2)) + abs(c(ilev,ind_4_4)) + eps) 1033 eps_4 = min(eps_4,1.0 _jprb)1081 eps_4 = min(eps_4,1.0) 1034 1082 1035 1083 mat(ind_4_2,ind_4_2) = mat(ind_4_2,ind_4_2) + indice_4(i4)%z1*v_4(ilev,i4)*(1. - eps_4)*c(ilev,ind_4_4)
Note: See TracChangeset
for help on using the changeset viewer.
