Changeset 7180 for branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90

Timestamp:

2016-11-03T16:41:10+01:00 (7 years ago)

Author:

aumont

Message:

various bug fixes on iron chemistry

File:

• branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90 (modified) (5 diffs)

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90

@@ -69 +69 @@
       INTEGER  ::   ji, jj, jk, jic, jn
       REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
-      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, fe3sol
+      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, fe3sol, fe3sol1
       REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag
       REAL(wp) ::   ztrc, zdust
@@ -279 +279 @@
                   zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
                ELSE
+                  zfeequi = zFe3(ji,jj,jk) * 1E-9
                   IF (ln_fecolloid) THEN
-                     zfeequi = zFe3(ji,jj,jk) * 1E-9
                      zhplus   = max( rtrn, hi(ji,jj,jk) )
-                     fe3sol  = fesol(ji,jj,jk,1) * ( fesol(ji,jj,jk,2) * zhplus**2  &
+                     fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
                      &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
                      &         + fesol(ji,jj,jk,5) / zhplus )
                      zfecoll = max( ( 0.1 * zFeL1(ji,jj,jk) * 1E-9 ), ( zFeL1(ji,jj,jk) * 1E-9 -fe3sol ) )
+#if defined key_ligand
+                     zligco  = max( ( 0.1 * trn(ji,jj,jk,jplgw) ), ( trn(ji,jj,jk,jplgw) - fe3sol ) )
+#endif 
                   ELSE
-                     zfeequi = zFe3(ji,jj,jk) * 1E-9 
                      zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
+#if defined key_ligand
+                     zligco  = 0.5 * trn(ji,jj,jk,jplgw)
+#endif 
                      fe3sol  = 0.
-                     kfep    = 0.
                   ENDIF
                ENDIF
@@ -324 +328 @@
                !  ----------------------------------------------------------------
                zlam1a  = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
-                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) + 0. * trb(ji,jj,jk,jppoc) )
+                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) )
                zaggdfea = zlam1a * zstep * zfecoll
-#if defined key_ligand
-               zligco = max( ( 0.1 * trn(ji,jj,jk,jplgw) ), ( trn(ji,jj,jk,jplgw) - fe3sol ) )
-               zaggliga = zlam1a * zstep * zligco
-#endif
                !
 #if defined key_kriest
                zaggdfeb = 0.
-               !
-#  if defined key_ligand
-               zaggligb = 0.
-#  endif
-               !
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea + zaggdfeb
-               !
 #else
-               !
                zlam1b = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
                zaggdfeb = zlam1b * zstep * zfecoll
-               !
-#  if defined key_ligand
-               zaggligb = zlam1b * zstep * zligco
-#  endif
+#endif
                ! precipitation of Fe3+, creation of nanoparticles
                precip(ji,jj,jk) = max( 0., (zfeequi - fe3sol) ) * kfep * zstep
@@ -355 +343 @@
                &                     - zcoag - precip(ji,jj,jk)
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
+#if ! defined key_kriest 
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
 #endif
 #if defined key_ligand
+               zaggliga = zlam1a * zstep * zligco
+#   if defined key_kriest
+               zaggligb = 0.
+#   else
+               zaggligb = zlam1b * zstep * zligco
+#   endif
                tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) + precip(ji,jj,jk)
                tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb
@@ -404 +399 @@
       ENDIF
       !
-                       CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
-      IF( ln_fechem )  CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
+      CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
+      IF( ln_fechem )  THEN
+         CALL wrk_dealloc( jpi, jpj,      zstrn, zstrn2 )
+         CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
+      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_fechem')