Changeset 1924 for trunk/LMDZ.TITAN/libf/phytitan

Timestamp:

May 2, 2018, 3:18:27 PM (7 years ago)

Author:

jvatant

Message:

Correct a if statement calling - safely - out-ouf-bounds arrays in calchim.F90.
Was a problem for exacting debug mode.
--JVO

File:

• trunk/LMDZ.TITAN/libf/phytitan/calchim.F90 (modified) (1 diff)

trunk/LMDZ.TITAN/libf/phytitan/calchim.F90

@@ -462 +462 @@
            ENDIF
 
-           DO i=1,nd_kim+1
-
-              krpddecm1 =   (   krpd(idec  ,i,ialt  ,klat)   * (1.0-factalt)                   &
-                              + krpd(idec  ,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
-                          * (   krpd(idec  ,i,ialt  ,klat+1) * (1.0-factalt)                   &
-                              + krpd(idec  ,i,ialt+1,klat+1) * factalt       ) * factlat
+           DO i=1,nd_kim+1 ! nd_kim+1 is dissociation of N2 by GCR
+
               krpddec   =   (   krpd(idec+1,i,ialt  ,klat)   * (1.0-factalt)                   &
                               + krpd(idec+1,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
                           * (   krpd(idec+1,i,ialt  ,klat+1) * (1.0-factalt)                   &
                               + krpd(idec+1,i,ialt+1,klat+1) * factalt       ) * factlat
-              krpddecp1 =   (   krpd(idec+2,i,ialt  ,klat)   * (1.0-factalt)                   &
-                              + krpd(idec+2,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
-                          * (   krpd(idec+2,i,ialt  ,klat+1) * (1.0-factalt)                   &
-                              + krpd(idec+2,i,ialt+1,klat+1) * factalt       ) * factlat
-
-              ! nd_kim+1 is dissociation of N2 by GCR
-              if ( factdec.lt.0. ) then 
+
+              if      ( factdec.lt.0. ) then 
+                 krpddecm1 =   (   krpd(idec  ,i,ialt  ,klat)   * (1.0-factalt)                   &
+                                 + krpd(idec  ,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
+                             * (   krpd(idec  ,i,ialt  ,klat+1) * (1.0-factalt)                   &
+                                 + krpd(idec  ,i,ialt+1,klat+1) * factalt       ) * factlat
                  krate(l,i) = krpddecm1 * abs(factdec) + krpddec   * ( 1.0 + factdec)
+              else if ( factdec.gt.0. ) then
+                 krpddecp1 =   (   krpd(idec+2,i,ialt  ,klat)   * (1.0-factalt)                   &
+                                 + krpd(idec+2,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
+                             * (   krpd(idec+2,i,ialt  ,klat+1) * (1.0-factalt)                   &
+                                 + krpd(idec+2,i,ialt+1,klat+1) * factalt       ) * factlat
+                 krate(l,i) = krpddecp1 * factdec      + krpddec   * ( 1.0 - factdec)
+              else if ( factdec.eq.0. ) then 
+                 krate(l,i) = krpddec
               endif
-              if ( factdec.gt.0. ) then
-                 krate(l,i) = krpddecp1 * factdec      + krpddec   * ( 1.0 - factdec)
-              endif
-              if ( factdec.eq.0. ) krate(l,i) = krpddec
 
            ENDDO ! i=1,nd_kim+1