source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/build/ppsrc/INCA_SRC/imp_slv.f90 @ 6610

Last change on this file since 6610 was 6610, checked in by acosce, 10 months ago

INCA used for ICOLMDZORINCA_CO2_Transport_GMD_2023
File size: 12.8 KB
Line 
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
65SUBROUTINE 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
80END SUBROUTINE IMP_SLV_INTI
81
82SUBROUTINE 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
338END SUBROUTINE IMP_SOL
339
Note: See TracBrowser for help on using the repository browser.