[6610] | 1 | |
---|
| 2 | |
---|
| 3 | |
---|
| 4 | |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | !$Id: imp_slv.F90 119 2009-03-27 10:06:33Z acosce $ |
---|
| 13 | !! ========================================================================= |
---|
| 14 | !! INCA - INteraction with Chemistry and Aerosols |
---|
| 15 | !! |
---|
| 16 | !! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE) |
---|
| 17 | !! Unite mixte CEA-CNRS-UVSQ |
---|
| 18 | !! |
---|
| 19 | !! Contributors to this INCA subroutine: |
---|
| 20 | !! |
---|
| 21 | !! Stacy Walters, NCAR, stacy@ucar.edu |
---|
| 22 | !! Didier Hauglustaine, LSCE. |
---|
| 23 | !! |
---|
| 24 | !! Anne Cozic, LSCE, anne.cozic@cea.fr |
---|
| 25 | !! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr |
---|
| 26 | !! |
---|
| 27 | !! This software is a computer program whose purpose is to simulate the |
---|
| 28 | !! atmospheric gas phase and aerosol composition. The model is designed to be |
---|
| 29 | !! used within a transport model or a general circulation model. This version |
---|
| 30 | !! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts |
---|
| 31 | !! for emissions, transport (resolved and sub-grid scale), photochemical |
---|
| 32 | !! transformations, and scavenging (dry deposition and washout) of chemical |
---|
| 33 | !! species and aerosols interactively in the GCM. Several versions of the INCA |
---|
| 34 | !! model are currently used depending on the envisaged applications with the |
---|
| 35 | !! chemistry-climate model. |
---|
| 36 | !! |
---|
| 37 | !! This software is governed by the CeCILL license under French law and |
---|
| 38 | !! abiding by the rules of distribution of free software. You can use, |
---|
| 39 | !! modify and/ or redistribute the software under the terms of the CeCILL |
---|
| 40 | !! license as circulated by CEA, CNRS and INRIA at the following URL |
---|
| 41 | !! "http://www.cecill.info". |
---|
| 42 | !! |
---|
| 43 | !! As a counterpart to the access to the source code and rights to copy, |
---|
| 44 | !! modify and redistribute granted by the license, users are provided only |
---|
| 45 | !! with a limited warranty and the software's author, the holder of the |
---|
| 46 | !! economic rights, and the successive licensors have only limited |
---|
| 47 | !! liability. |
---|
| 48 | !! |
---|
| 49 | !! In this respect, the user's attention is drawn to the risks associated |
---|
| 50 | !! with loading, using, modifying and/or developing or reproducing the |
---|
| 51 | !! software by the user in light of its specific status of free software, |
---|
| 52 | !! that may mean that it is complicated to manipulate, and that also |
---|
| 53 | !! therefore means that it is reserved for developers and experienced |
---|
| 54 | !! professionals having in-depth computer knowledge. Users are therefore |
---|
| 55 | !! encouraged to load and test the software's suitability as regards their |
---|
| 56 | !! requirements in conditions enabling the security of their systems and/or |
---|
| 57 | !! data to be ensured and, more generally, to use and operate it in the |
---|
| 58 | !! same conditions as regards security. |
---|
| 59 | !! |
---|
| 60 | !! The fact that you are presently reading this means that you have had |
---|
| 61 | !! knowledge of the CeCILL license and that you accept its terms. |
---|
| 62 | !! ========================================================================= |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | SUBROUTINE IMP_SLV_INTI() |
---|
| 66 | !----------------------------------------------------------------------- |
---|
| 67 | ! ... Initialize the implict solver |
---|
| 68 | !----------------------------------------------------------------------- |
---|
| 69 | |
---|
| 70 | USE IMP_SLV0 |
---|
| 71 | |
---|
| 72 | IMPLICIT NONE |
---|
| 73 | |
---|
| 74 | |
---|
| 75 | factor(:) = .true. |
---|
| 76 | |
---|
| 77 | EPSILON(:) = .001 |
---|
| 78 | EPSILON((/1,3,4,7,8,9,10,11,12/)) = .0001 |
---|
| 79 | |
---|
| 80 | END SUBROUTINE IMP_SLV_INTI |
---|
| 81 | |
---|
| 82 | SUBROUTINE IMP_SOL( base_sol & |
---|
| 83 | , reaction_rates & |
---|
| 84 | , het_rates & |
---|
| 85 | , extfrc & |
---|
| 86 | , o3_prod & |
---|
| 87 | , o3_loss & |
---|
| 88 | , pdel & |
---|
| 89 | , hnm & |
---|
| 90 | , nstep & |
---|
| 91 | , delt_inverse & |
---|
| 92 | , lat , wetloss) |
---|
| 93 | !----------------------------------------------------------------------- |
---|
| 94 | ! ... Imp_sol advances the volumetric mixing ratio |
---|
| 95 | ! forward one time step via the fully implicit |
---|
| 96 | ! Euler scheme |
---|
| 97 | ! Stacy Walters, NCAR, 1999. |
---|
| 98 | !----------------------------------------------------------------------- |
---|
| 99 | USE CHEM_MODS, ONLY: clsmap, permute, diag_map, adv_mass |
---|
| 100 | USE SPECIES_NAMES |
---|
| 101 | USE IMP_SLV0 |
---|
| 102 | USE CHEM_TRACNM |
---|
| 103 | USE INCA_DIM |
---|
| 104 | USE RATE_INDEX_MOD |
---|
| 105 | |
---|
| 106 | IMPLICIT NONE |
---|
| 107 | |
---|
| 108 | !----------------------------------------------------------------------- |
---|
| 109 | ! ... Dummy args |
---|
| 110 | !----------------------------------------------------------------------- |
---|
| 111 | INTEGER, INTENT(in) :: nstep ! time step index (zero based) |
---|
| 112 | INTEGER, INTENT(in) :: lat ! lat index (1...#PLAT) |
---|
| 113 | REAL, INTENT(in) :: delt_inverse ! 1. / delt ( delt == time step in seconds) |
---|
| 114 | REAL, INTENT(in) :: reaction_rates(PLNPLV,12) |
---|
| 115 | REAL, INTENT(inout) :: o3_prod(PLNPLV), o3_loss(PLNPLV) |
---|
| 116 | REAL, INTENT(inout) :: hnm(PLNPLV) |
---|
| 117 | REAL, INTENT(in) :: pdel(PLNPLV) ! delta press across midpoints |
---|
| 118 | REAL, INTENT(in) :: het_rates(PLNPLV,1) |
---|
| 119 | REAL, INTENT(in) :: extfrc(PLNPLV,1) |
---|
| 120 | REAL, INTENT(inout) :: base_sol(PLNPLV,8) |
---|
| 121 | REAL, INTENT(out) :: wetloss(PLNPLV,1) |
---|
| 122 | !----------------------------------------------------------------------- |
---|
| 123 | ! ... Local variables |
---|
| 124 | !----------------------------------------------------------------------- |
---|
| 125 | INTEGER :: nr_iter |
---|
| 126 | INTEGER :: lev |
---|
| 127 | INTEGER :: ofl, ofu |
---|
| 128 | INTEGER :: i,j,k,m,n |
---|
| 129 | INTEGER :: hndx |
---|
| 130 | REAL :: max_delta(0) |
---|
| 131 | REAL, DIMENSION(PLON,0) :: sys_jac |
---|
| 132 | REAL, DIMENSION(PLON,0) :: lin_jac |
---|
| 133 | REAL, DIMENSION(PLON,0) ::solution |
---|
| 134 | REAL, DIMENSION(PLON,0) ::forcing |
---|
| 135 | REAL, DIMENSION(PLON,0) ::iter_invariant |
---|
| 136 | REAL, DIMENSION(PLON,0) ::prod |
---|
| 137 | REAL, DIMENSION(PLON,0) ::loss |
---|
| 138 | REAL, DIMENSION(PLON) :: wrk |
---|
| 139 | REAL, PARAMETER :: undef = 1.e-40 |
---|
| 140 | REAL, DIMENSION(PLNPLV,0) :: ind_prd |
---|
| 141 | REAL :: fac ! factor for computation of wet dep of gases |
---|
| 142 | LOGICAL :: converged(0) |
---|
| 143 | LOGICAL :: converged_all(PLON,0) |
---|
| 144 | REAL :: saved_base_sol(PLON,8) |
---|
| 145 | !----------------------------------------------------------------------- |
---|
| 146 | ! ... Function declarations |
---|
| 147 | !----------------------------------------------------------------------- |
---|
| 148 | |
---|
| 149 | max_delta(:) = 0. |
---|
| 150 | |
---|
| 151 | !----------------------------------------------------------------------- |
---|
| 152 | ! ... Class independent forcing |
---|
| 153 | !----------------------------------------------------------------------- |
---|
| 154 | CALL INDPRD( 4, & |
---|
| 155 | ind_prd, & |
---|
| 156 | base_sol, & |
---|
| 157 | extfrc, & |
---|
| 158 | reaction_rates ) |
---|
| 159 | |
---|
| 160 | DO lev = 1,PLEV |
---|
| 161 | saved_base_sol(:,:) = 0. |
---|
| 162 | converged_all(:,:) = .FALSE. |
---|
| 163 | converged(:) = .FALSE. |
---|
| 164 | ofl = (lev - 1)*PLON + 1 |
---|
| 165 | ofu = ofl + PLON - 1 |
---|
| 166 | !----------------------------------------------------------------------- |
---|
| 167 | ! ... Transfer from base to class array |
---|
| 168 | !----------------------------------------------------------------------- |
---|
| 169 | DO k = 1,0 |
---|
| 170 | j = clsmap(k,4) |
---|
| 171 | m = permute(k,4) |
---|
| 172 | solution(:,m) = base_sol(ofl:ofu,j) |
---|
| 173 | END DO ! enddo k |
---|
| 174 | |
---|
| 175 | !----------------------------------------------------------------------- |
---|
| 176 | ! ... Set the iteration invariant part of the function F(y) |
---|
| 177 | ! ... If there is "independent" production put it in the forcing |
---|
| 178 | !----------------------------------------------------------------------- |
---|
| 179 | DO m = 1,0 |
---|
| 180 | iter_invariant(:,m) = delt_inverse * solution(:,m) & |
---|
| 181 | + ind_prd(ofl:ofu,m) |
---|
| 182 | END DO ! enddo m |
---|
| 183 | |
---|
| 184 | !----------------------------------------------------------------------- |
---|
| 185 | ! ... The linear component |
---|
| 186 | !----------------------------------------------------------------------- |
---|
| 187 | DO m = 1,0 |
---|
| 188 | DO k = 1,PLON |
---|
| 189 | lin_jac(k,m) = 0. |
---|
| 190 | END DO ! enddo k |
---|
| 191 | END DO ! enddo m |
---|
| 192 | DO j = 1,0 |
---|
| 193 | m = diag_map(j) |
---|
| 194 | lin_jac(:,m) = -delt_inverse |
---|
| 195 | END DO ! END DO j |
---|
| 196 | |
---|
| 197 | !======================================================================= |
---|
| 198 | ! The Newton-Raphson iteration for F(y) = 0 |
---|
| 199 | !======================================================================= |
---|
| 200 | DO nr_iter = 1,itermax |
---|
| 201 | !----------------------------------------------------------------------- |
---|
| 202 | ! ... The non-linear component |
---|
| 203 | !----------------------------------------------------------------------- |
---|
| 204 | IF( factor(nr_iter) ) THEN |
---|
| 205 | DO m = 1,0 |
---|
| 206 | DO k = 1,PLON |
---|
| 207 | sys_jac(k,m) = lin_jac(k,m) |
---|
| 208 | END DO ! END DO k |
---|
| 209 | END DO ! END DO m |
---|
| 210 | !----------------------------------------------------------------------- |
---|
| 211 | ! ... Factor the "system" matrix |
---|
| 212 | !----------------------------------------------------------------------- |
---|
| 213 | CALL LU_FAC( sys_jac ) |
---|
| 214 | END IF |
---|
| 215 | !----------------------------------------------------------------------- |
---|
| 216 | ! ... Form F(y) |
---|
| 217 | !----------------------------------------------------------------------- |
---|
| 218 | CALL IMP_PROD_LOSS( prod, & |
---|
| 219 | loss, & |
---|
| 220 | base_sol(ofl,1), & |
---|
| 221 | reaction_rates(ofl,1), & |
---|
| 222 | het_rates(ofl,1) ) |
---|
| 223 | |
---|
| 224 | |
---|
| 225 | DO m = 1,0 |
---|
| 226 | DO k = 1,PLON |
---|
| 227 | forcing(k,m) = solution(k,m)*delt_inverse & |
---|
| 228 | - (iter_invariant(k,m) + prod(k,m) - loss(k,m)) |
---|
| 229 | |
---|
| 230 | END DO ! END DO k |
---|
| 231 | END DO ! END DO m |
---|
| 232 | |
---|
| 233 | |
---|
| 234 | !----------------------------------------------------------------------- |
---|
| 235 | ! ... Solve for the mixing ratio at t(n+1) |
---|
| 236 | !----------------------------------------------------------------------- |
---|
| 237 | CALL LU_SLV( sys_jac, forcing ) |
---|
| 238 | DO m = 1,0 |
---|
| 239 | DO k = 1,PLON |
---|
| 240 | solution(k,m) = solution(k,m) + forcing(k,m) |
---|
| 241 | END DO ! END DO k |
---|
| 242 | END DO ! END do m |
---|
| 243 | |
---|
| 244 | !----------------------------------------------------------------------- |
---|
| 245 | ! ... Convergence measures |
---|
| 246 | !----------------------------------------------------------------------- |
---|
| 247 | IF( nr_iter > 1 ) THEN |
---|
| 248 | DO k = 1,0 |
---|
| 249 | m = permute(k,4) |
---|
| 250 | WHERE( ABS(solution(:,m)) > undef ) |
---|
| 251 | wrk(:) = ABS( forcing(:,m)/solution(:,m) ) |
---|
| 252 | ELSEWHERE |
---|
| 253 | wrk(:) = undef |
---|
| 254 | END WHERE |
---|
| 255 | max_delta(k) = MAXVAL( wrk(:) ) |
---|
| 256 | END DO ! END DO k |
---|
| 257 | END IF |
---|
| 258 | |
---|
| 259 | !----------------------------------------------------------------------- |
---|
| 260 | ! ... Limit iterate |
---|
| 261 | !----------------------------------------------------------------------- |
---|
| 262 | DO m = 1,0 |
---|
| 263 | WHERE( solution(:,m) < 0. ) |
---|
| 264 | solution(:,m) = 0. |
---|
| 265 | END WHERE |
---|
| 266 | END DO ! END DO m |
---|
| 267 | |
---|
| 268 | !----------------------------------------------------------------------- |
---|
| 269 | ! ... Transfer latest solution back to "base" array |
---|
| 270 | !----------------------------------------------------------------------- |
---|
| 271 | DO k = 1,0 |
---|
| 272 | j = clsmap(k,4) |
---|
| 273 | m = permute(k,4) |
---|
| 274 | base_sol(ofl:ofu,j) = solution(:,m) |
---|
| 275 | END DO !END DO k |
---|
| 276 | |
---|
| 277 | !----------------------------------------------------------------------- |
---|
| 278 | ! Check for convergence |
---|
| 279 | !----------------------------------------------------------------------- |
---|
| 280 | IF( nr_iter > 1 ) THEN |
---|
| 281 | DO k = 1,0 |
---|
| 282 | m = permute(k,4) |
---|
| 283 | converged(k) = .NOT. ANY( ABS(forcing(:,m))& |
---|
| 284 | > EPSILON(k)*ABS(solution(:,m)) ) |
---|
| 285 | |
---|
| 286 | j = clsmap(k,4) |
---|
| 287 | DO i=1,PLON |
---|
| 288 | IF ( (.NOT. converged_all(i,k)) .AND. & |
---|
| 289 | (ABS(forcing(i,m)) < EPSILON(k)*ABS(solution(i,m)))) THEN |
---|
| 290 | saved_base_sol(i,j)=solution(i,m) |
---|
| 291 | converged_all(i,k) = .TRUE. |
---|
| 292 | ENDIF |
---|
| 293 | ENDDO |
---|
| 294 | |
---|
| 295 | END DO ! END DO k |
---|
| 296 | IF( ALL( converged(:) ) ) THEN |
---|
| 297 | EXIT |
---|
| 298 | END IF |
---|
| 299 | END IF |
---|
| 300 | END DO ! END DO nr_iter |
---|
| 301 | |
---|
| 302 | IF( ANY( .NOT. converged(:) ) ) THEN |
---|
| 303 | |
---|
| 304 | ! WRITE(*,*) 'IMP_SOL: Failed to converge @ (lev) = ',lev |
---|
| 305 | DO m = 1,0 |
---|
| 306 | ! IF( .NOT. converged(m) ) THEN |
---|
| 307 | ! WRITE(*,'(1x,a8,1x,1pe10.3)') solsym(clsmap(m,4)), max_delta(m) |
---|
| 308 | ! ENDIF |
---|
| 309 | DO i=1,PLON |
---|
| 310 | IF (converged_all(i,m)) THEN |
---|
| 311 | j = clsmap(m,4) |
---|
| 312 | base_sol(ofl+i-1,j) = saved_base_sol(i,j) |
---|
| 313 | ENDIF |
---|
| 314 | ENDDO |
---|
| 315 | END DO !END DO m |
---|
| 316 | |
---|
| 317 | ELSE |
---|
| 318 | DO k = 1,0 |
---|
| 319 | j = clsmap(k,4) |
---|
| 320 | base_sol(ofl:ofu,j) = saved_base_sol(:,j) |
---|
| 321 | END DO !END DO k |
---|
| 322 | END IF |
---|
| 323 | |
---|
| 324 | |
---|
| 325 | !----------------------------------------------------------------------- |
---|
| 326 | ! ... Ozone prod and loss terms |
---|
| 327 | !----------------------------------------------------------------------- |
---|
| 328 | DO k = ofl,ofu |
---|
| 329 | !----------------------------------------------------------------------- |
---|
| 330 | ! ... Ox production (only valid for the troposphere!) |
---|
| 331 | !----------------------------------------------------------------------- |
---|
| 332 | |
---|
| 333 | |
---|
| 334 | END DO ! END DO k |
---|
| 335 | |
---|
| 336 | |
---|
| 337 | END DO ! END DO lev |
---|
| 338 | END SUBROUTINE IMP_SOL |
---|
| 339 | |
---|