Changeset 6453 for branches/CNRS/dev_r6270_PISCES_QUOTA

Timestamp:

2016-04-08T10:10:11+02:00 (8 years ago)

Author:

aumont

Message:

New developments of PISCES (QUOTA, ligands, lability, ...)

Location:

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM

Files:

• CONFIG/ORCA2_OFF_PISCES_QUOTA (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00 (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/domain_def.xml (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/field_def.xml (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/iodef.xml (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/namelist_cfg (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/namelist_pisces_cfg (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/namelist_pisces_ref (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/namelist_ref (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/namelist_top_cfg (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/EXP00/namelist_top_ref (added)
• CONFIG/ORCA2_OFF_PISCES_QUOTA/cpp_ORCA2_OFF_PISCES_QUOTA.fcm (added)
• NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90 (added)
• NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90 (modified) (2 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90 (modified) (13 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90 (modified) (17 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90 (modified) (19 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90 (modified) (2 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90 (added)
• NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90 (modified) (3 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90 (modified) (10 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90 (modified) (4 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90 (modified) (3 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90 (modified) (3 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90 (modified) (13 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90 (added)
• NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90 (modified) (3 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90 (modified) (12 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90 (modified) (16 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90 (modified) (9 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90 (modified) (16 diffs)
• NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90 (modified) (9 diffs)
• NEMO/TOP_SRC/PISCES/P5Z (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zagg.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zbio.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zlim.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zmeso.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zmicro.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zmort.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zrem.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zsed.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zsink.F90 (added)
• NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90 (added)
• NEMO/TOP_SRC/PISCES/par_pisces.F90 (modified) (5 diffs)
• NEMO/TOP_SRC/PISCES/sms_pisces.F90 (modified) (8 diffs)
• NEMO/TOP_SRC/PISCES/trcice_pisces.F90 (modified) (5 diffs)
• NEMO/TOP_SRC/PISCES/trcini_pisces.F90 (modified) (11 diffs)
• NEMO/TOP_SRC/PISCES/trcnam_pisces.F90 (modified) (3 diffs)
• NEMO/TOP_SRC/PISCES/trcsms_pisces.F90 (modified) (2 diffs)
• NEMO/TOP_SRC/PISCES/trcwri_pisces.F90 (modified) (3 diffs)

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

@@ -24 +24 @@
    USE p4zmicro        !  Sources and sinks of microzooplankton
    USE p4zmeso         !  Sources and sinks of mesozooplankton
-   USE p4zrem          !  Remineralisation of organic matter
-   USE p4zfechem
+   USE p4zrem          !  Remineralisation of organic/inorganic matter
+   USE p4zpoc          !  Remineralization of organic particles
+   USE p4zagg          !  Aggregation of particles
+   USE p4zlys          !  Dissolution of calcite
+   USE p4zfechem       !  Iron chemistry
+   USE p4zligand
    USE prtctl_trc      !  print control for debugging
    USE iom             !  I/O manager
@@ -76 +80 @@
 
           
-      CALL p4z_opt  ( kt, knt )     ! Optic: PAR in the water column
-      CALL p4z_sink ( kt, knt )     ! vertical flux of particulate organic matter
-      CALL p4z_fechem(kt, knt )     ! Iron chemistry/scavenging
-      CALL p4z_lim  ( kt, knt )     ! co-limitations by the various nutrients
-      CALL p4z_prod ( kt, knt )     ! phytoplankton growth rate over the global ocean. 
-      !                             ! (for each element : C, Si, Fe, Chl )
-      CALL p4z_mort ( kt      )     ! phytoplankton mortality
-     !                             ! zooplankton sources/sinks routines 
-      CALL p4z_micro( kt, knt )           ! microzooplankton
-      CALL p4z_meso ( kt, knt )           ! mesozooplankton
-      CALL p4z_rem  ( kt, knt )     ! remineralization terms of organic matter+scavenging of Fe
+      CALL p4z_opt   ( kt, knt )     ! Optic: PAR in the water column
+      CALL p4z_sink  ( kt, knt )     ! vertical flux of particulate organic matter
+      CALL p4z_lys   (kt, knt )      ! Dissolution of calcite
+      CALL p4z_fechem(kt, knt )      ! Iron chemistry/scavenging
+      CALL p4z_lim   ( kt, knt )     ! co-limitations by the various nutrients
+      CALL p4z_prod  ( kt, knt )     ! phytoplankton growth rate over the global ocean. 
+      !                              ! (for each element : C, Si, Fe, Chl )
+      CALL p4z_mort  ( kt      )     ! phytoplankton mortality
+      !                              ! zooplankton sources/sinks routines 
+      CALL p4z_micro ( kt, knt )           ! microzooplankton
+      CALL p4z_meso  ( kt, knt )           ! mesozooplankton
+      CALL p4z_agg   ( kt, knt )     ! Aggregation of particles
+      CALL p4z_rem   ( kt, knt )     ! remineralization terms of organic matter+scavenging of Fe
+      CALL p4z_poc   ( kt, knt )     ! Remineralization of organic particles
+      CALL p4z_ligand(kt, knt )
       !                             ! test if tracers concentrations fall below 0.
       !                                                             !

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

@@ -11 +11 @@
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
    !!                  !  2011-02  (J. Simeon, J.Orr ) update O2 solubility constants
+   !!             3.6  !  2016-03  (O. Aumont) Change chemistry to MOCSY standards
    !!----------------------------------------------------------------------
-#if defined key_pisces
+#if defined key_pisces || defined key_pisces_quota
    !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
+   !!   'key_pisces*'                                      PISCES bio-model
    !!----------------------------------------------------------------------
    !!   p4z_che      :  Sea water chemistry computed following OCMIP protocol
@@ -26 +27 @@
    PRIVATE
 
-   PUBLIC   p4z_che         !
-   PUBLIC   p4z_che_alloc   !
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sio3eq   ! chemistry of Si
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fekeq    ! chemistry of Fe
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemc    ! Solubilities of O2 and CO2
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2    ! Solubilities of O2 and CO2
+   PUBLIC   p4z_che          !
+   PUBLIC   p4z_che_alloc    !
+   PUBLIC   ahini_for_at     !
+   PUBLIC   solve_at_general !
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: sio3eq   ! chemistry of Si
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: fekeq    ! chemistry of Fe
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: chemc    ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: chemo2    ! Solubilities of O2 and CO2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: fesol    ! solubility of Fe
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akw3       !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akf3       !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aks3       !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak1p3      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak2p3      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak3p3      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksi3      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   borat      !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fluorid    !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sulfat     !: ???
+
+   !!* Variable for chemistry of the CO2 cycle
 
    REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm
 
-   REAL(wp) ::   salchl = 1. / 1.80655    ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)
    REAL(wp) ::   o2atm  = 1. / ( 1000. * 0.20946 )  
 
-   REAL(wp) ::   akcc1  = -171.9065       ! coeff. for apparent solubility equilibrium
-   REAL(wp) ::   akcc2  =   -0.077993     ! Millero et al. 1995 from Mucci 1983
-   REAL(wp) ::   akcc3  = 2839.319        
-   REAL(wp) ::   akcc4  =   71.595        
-   REAL(wp) ::   akcc5  =   -0.77712      
-   REAL(wp) ::   akcc6  =    0.00284263   
-   REAL(wp) ::   akcc7  =  178.34        
-   REAL(wp) ::   akcc8  =   -0.07711     
-   REAL(wp) ::   akcc9  =    0.0041249   
-
-   REAL(wp) ::   rgas   = 83.143         ! universal gas constants
-   REAL(wp) ::   oxyco  = 1. / 22.4144   ! converts from liters of an ideal gas to moles
-
-   REAL(wp) ::   bor1   = 0.00023        ! borat constants
-   REAL(wp) ::   bor2   = 1. / 10.82
-
-   REAL(wp) ::   ca0    = -162.8301      ! WEISS & PRICE 1980, units mol/(kg atm)
-   REAL(wp) ::   ca1    =  218.2968
-   REAL(wp) ::   ca2    =   90.9241
-   REAL(wp) ::   ca3    =   -1.47696
-   REAL(wp) ::   ca4    =    0.025695
-   REAL(wp) ::   ca5    =   -0.025225
-   REAL(wp) ::   ca6    =    0.0049867
-
-   REAL(wp) ::   c10    = -3670.7        ! Coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)   
-   REAL(wp) ::   c11    =    62.008     
-   REAL(wp) ::   c12    =    -9.7944    
-   REAL(wp) ::   c13    =     0.0118     
-   REAL(wp) ::   c14    =    -0.000116
-
-   REAL(wp) ::   c20    = -1394.7       ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)   
-   REAL(wp) ::   c21    =    -4.777   
-   REAL(wp) ::   c22    =     0.0184   
-   REAL(wp) ::   c23    =    -0.000118
-
-   REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate
-   REAL(wp) ::   st2    =  1./96.062    !  (Morris & Riley 1966)
-   REAL(wp) ::   ks0    =    141.328 
-   REAL(wp) ::   ks1    =  -4276.1  
-   REAL(wp) ::   ks2    =    -23.093
-   REAL(wp) ::   ks3    = -13856.  
-   REAL(wp) ::   ks4    =   324.57 
-   REAL(wp) ::   ks5    =   -47.986
-   REAL(wp) ::   ks6    =  35474. 
-   REAL(wp) ::   ks7    =   -771.54
-   REAL(wp) ::   ks8    =    114.723
-   REAL(wp) ::   ks9    =  -2698.  
-   REAL(wp) ::   ks10   =   1776. 
-   REAL(wp) ::   ks11   =      1.
-   REAL(wp) ::   ks12   =     -0.001005 
-
-   REAL(wp) ::   ft1    =    0.000067   ! constants for calculate concentrations for fluorides
-   REAL(wp) ::   ft2    = 1./18.9984    ! (Dickson & Riley 1979 )
-   REAL(wp) ::   kf0    =  -12.641    
-   REAL(wp) ::   kf1    = 1590.2    
-   REAL(wp) ::   kf2    =    1.525    
-   REAL(wp) ::   kf3    =    1.0     
-   REAL(wp) ::   kf4    =   -0.001005
-
-   REAL(wp) ::   cb0    = -8966.90      ! Coeff. for 1. dissoc. of boric acid 
-   REAL(wp) ::   cb1    = -2890.53      ! (Dickson and Goyet, 1994)
-   REAL(wp) ::   cb2    =   -77.942
-   REAL(wp) ::   cb3    =     1.728
-   REAL(wp) ::   cb4    =    -0.0996
-   REAL(wp) ::   cb5    =   148.0248
-   REAL(wp) ::   cb6    =   137.1942
-   REAL(wp) ::   cb7    =     1.62142
-   REAL(wp) ::   cb8    =   -24.4344
-   REAL(wp) ::   cb9    =   -25.085
-   REAL(wp) ::   cb10   =    -0.2474 
-   REAL(wp) ::   cb11   =     0.053105
-
-   REAL(wp) ::   cw0    = -13847.26     ! Coeff. for dissoc. of water (Dickson and Riley, 1979 )
-   REAL(wp) ::   cw1    =    148.9652  
-   REAL(wp) ::   cw2    =    -23.6521
-   REAL(wp) ::   cw3    =    118.67 
-   REAL(wp) ::   cw4    =     -5.977 
-   REAL(wp) ::   cw5    =      1.0495  
-   REAL(wp) ::   cw6    =     -0.01615
-
-   !                                    ! volumetric solubility constants for o2 in ml/L  
-   REAL(wp) ::   ox0    =  2.00856      ! from Table 1 for Eq 8 of Garcia and Gordon, 1992.
-   REAL(wp) ::   ox1    =  3.22400      ! corrects for moisture and fugacity, but not total atmospheric pressure
-   REAL(wp) ::   ox2    =  3.99063      !      Original PISCES code noted this was a solubility, but 
-   REAL(wp) ::   ox3    =  4.80299      ! was in fact a bunsen coefficient with units L-O2/(Lsw atm-O2)
-   REAL(wp) ::   ox4    =  9.78188e-1   ! Hence, need to divide EXP( zoxy ) by 1000, ml-O2 => L-O2
-   REAL(wp) ::   ox5    =  1.71069      ! and atcox = 0.20946 to add the 1/atm dimension.
-   REAL(wp) ::   ox6    = -6.24097e-3   
-   REAL(wp) ::   ox7    = -6.93498e-3 
-   REAL(wp) ::   ox8    = -6.90358e-3
-   REAL(wp) ::   ox9    = -4.29155e-3 
-   REAL(wp) ::   ox10   = -3.11680e-7 
+   REAL(wp) ::   rgas   = 83.14472      ! universal gas constants
+   REAL(wp) ::   oxyco  = 1. / 22.4144  ! converts from liters of an ideal gas to moles
 
    !                                    ! coeff. for seawater pressure correction : millero 95
    !                                    ! AGRIF doesn't like the DATA instruction
-   REAL(wp) :: devk11  = -25.5
-   REAL(wp) :: devk12  = -15.82
-   REAL(wp) :: devk13  = -29.48
-   REAL(wp) :: devk14  = -25.60
-   REAL(wp) :: devk15  = -48.76
+   REAL(wp) :: devk10  = -25.5
+   REAL(wp) :: devk11  = -15.82
+   REAL(wp) :: devk12  = -29.48
+   REAL(wp) :: devk13  = -20.02
+   REAL(wp) :: devk14  = -18.03
+   REAL(wp) :: devk15  = -9.78
+   REAL(wp) :: devk16  = -48.76
+   REAL(wp) :: devk17  = -14.51
+   REAL(wp) :: devk18  = -23.12
+   REAL(wp) :: devk19  = -26.57
+   REAL(wp) :: devk110  = -29.48
    !
-   REAL(wp) :: devk21  = 0.1271
-   REAL(wp) :: devk22  = -0.0219
-   REAL(wp) :: devk23  = 0.1622
-   REAL(wp) :: devk24  = 0.2324
-   REAL(wp) :: devk25  = 0.5304
+   REAL(wp) :: devk20  = 0.1271
+   REAL(wp) :: devk21  = -0.0219
+   REAL(wp) :: devk22  = 0.1622
+   REAL(wp) :: devk23  = 0.1119
+   REAL(wp) :: devk24  = 0.0466
+   REAL(wp) :: devk25  = -0.0090
+   REAL(wp) :: devk26  = 0.5304
+   REAL(wp) :: devk27  = 0.1211
+   REAL(wp) :: devk28  = 0.1758
+   REAL(wp) :: devk29  = 0.2020
+   REAL(wp) :: devk210  = 0.1622
    !
+   REAL(wp) :: devk30  = 0.
    REAL(wp) :: devk31  = 0.
-   REAL(wp) :: devk32  = 0.
-   REAL(wp) :: devk33  = 2.608E-3
-   REAL(wp) :: devk34  = -3.6246E-3
-   REAL(wp) :: devk35  = 0.
+   REAL(wp) :: devk32  = 2.608E-3
+   REAL(wp) :: devk33  = -1.409e-3
+   REAL(wp) :: devk34  = 0.316e-3
+   REAL(wp) :: devk35  = -0.942e-3
+   REAL(wp) :: devk36  = 0.
+   REAL(wp) :: devk37  = -0.321e-3
+   REAL(wp) :: devk38  = -2.647e-3
+   REAL(wp) :: devk39  = -3.042e-3
+   REAL(wp) :: devk310  = -2.6080e-3
    !
-   REAL(wp) :: devk41  = -3.08E-3
-   REAL(wp) :: devk42  = 1.13E-3
-   REAL(wp) :: devk43  = -2.84E-3
-   REAL(wp) :: devk44  = -5.13E-3
-   REAL(wp) :: devk45  = -11.76E-3
+   REAL(wp) :: devk40  = -3.08E-3
+   REAL(wp) :: devk41  = 1.13E-3
+   REAL(wp) :: devk42  = -2.84E-3
+   REAL(wp) :: devk43  = -5.13E-3
+   REAL(wp) :: devk44  = -4.53e-3
+   REAL(wp) :: devk45  = -3.91e-3
+   REAL(wp) :: devk46  = -11.76e-3
+   REAL(wp) :: devk47  = -2.67e-3
+   REAL(wp) :: devk48  = -5.15e-3
+   REAL(wp) :: devk49  = -4.08e-3
+   REAL(wp) :: devk410  = -2.84e-3
    !
-   REAL(wp) :: devk51  = 0.0877E-3
-   REAL(wp) :: devk52  = -0.1475E-3     
-   REAL(wp) :: devk53  = 0.
-   REAL(wp) :: devk54  = 0.0794E-3      
-   REAL(wp) :: devk55  = 0.3692E-3      
+   REAL(wp) :: devk50  = 0.0877E-3
+   REAL(wp) :: devk51  = -0.1475E-3     
+   REAL(wp) :: devk52  = 0.
+   REAL(wp) :: devk53  = 0.0794E-3      
+   REAL(wp) :: devk54  = 0.09e-3
+   REAL(wp) :: devk55  = 0.054e-3
+   REAL(wp) :: devk56  = 0.3692E-3
+   REAL(wp) :: devk57  = 0.0427e-3
+   REAL(wp) :: devk58  = 0.09e-3
+   REAL(wp) :: devk59  = 0.0714e-3
+   REAL(wp) :: devk510  = 0.0
+   !
+
+
+   ! General parameters
+   REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp
+   REAL(wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp
+
+   ! Maximum number of iterations for each method
+   INTEGER, PARAMETER :: jp_maxniter_atgen    = 20
+
+   ! Bookkeeping variables for each method
+   ! - SOLVE_AT_GENERAL
+   INTEGER :: niter_atgen    = jp_maxniter_atgen
+
+
 
    !!* Substitution
@@ -182 +154 @@
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   ztkel, zt   , zt2   , zsal  , zsal2 , zbuf1 , zbuf2
+      REAL(wp) ::   ztkel, ztkel1, zt , zt2   , zsal  , zsal2 , zbuf1 , zbuf2
       REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
       REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
-      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1
+      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat
       REAL(wp) ::   zis  , zis2 , zsal15, zisqrt
       REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
+      REAL(wp) ::   zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi
       REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1
+      REAL(wp) ::   total2free, free2SWS, total2SWS, SWS2total
+
       !!---------------------------------------------------------------------
       !
@@ -200 +175 @@
          DO ji = 1, jpi
             !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tsn(ji,jj,1,jp_tem) + 273.16
+            ztkel = tsn(ji,jj,1,jp_tem) + 273.15
             zt    = ztkel * 0.01
-            zt2   = zt * zt
             zsal  = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35.
-            zsal2 = zsal * zsal
-            zlogt = LOG( zt )
             !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
             !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
-            zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 )
-            !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 in ml/L (EQ. 8, GARCIA AND GORDON, 1992)
-            ztgg  = LOG( ( 298.15 - tsn(ji,jj,1,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
-            ztgg2 = ztgg  * ztgg
-            ztgg3 = ztgg2 * ztgg
-            ztgg4 = ztgg3 * ztgg
-            ztgg5 = ztgg4 * ztgg
-            zoxy  = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5   &
-                   + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) +  ox10 * zsal2
-
+            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    &
+            &       + 0.0047036e-4*ztkel**2)
             !                             ! SET SOLUBILITIES OF O2 AND CO2 
-            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm)
-            chemc(ji,jj,2) = ( EXP( zoxy  ) * o2atm ) * oxyco              ! mol/(L atm)
+            chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 ! mol/(kg atm)
+            chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3
+            chemc(ji,jj,3) = 57.7 - 0.118*ztkel
             !
          END DO
@@ -233 +198 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-              ztkel = tsn(ji,jj,jk,jp_tem) + 273.16
+              ztkel = tsn(ji,jj,jk,jp_tem) + 273.15
               zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35.
               zsal2 = zsal * zsal
@@ -241 +206 @@
               ztgg4 = ztgg3 * ztgg
               ztgg5 = ztgg4 * ztgg
-              zoxy  = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5   &
-                     + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) +  ox10 * zsal2
+
+              zoxy  = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3    &
+              &       + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3   &
+              &       - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 )   &
+              &       - 3.11680e-7 * zsal2
               chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox     ! mol/(L atm)
             END DO
@@ -259 +227 @@
             DO ji = 1, jpi
 
-               ! SET PRESSION
-               zpres   = 1.025e-1 * fsdept(ji,jj,jk)
+               ! SET PRESSION ACCORDING TO SAUNDER (1980)
+               zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
+               zc1 = 5.92E-3 + zplat**2 * 5.25E-3
+               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*fsdept(ji,jj,jk)))) / 4.42E-6
+               zpres = zpres / 10.0
 
                ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.16
+               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.15
                zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
@@ -275 +246 @@
 
                ! CHLORINITY (WOOSTER ET AL., 1969)
-               zcl     = zsal * salchl
+               zcl     = zsal / 1.80655
 
                ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
-               zst     = st1 * zcl * st2
+               zst     = 0.14 * zcl /96.062
 
                ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
-               zft     = ft1 * zcl * ft2
+               zft     = 0.000067 * zcl /18.9984
 
                ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
-               zcks    = EXP(  ks1 * ztr + ks0 + ks2 * zlogt                           &
-                  &                     + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt   &
-                  &                     + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis      &
-                  &                     + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal )  )
+               zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         &
+               &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt &
+               &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    &
+               &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         &
+               &         + LOG(1.0 - 0.001005 * zsal))
+
 
                ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
-               zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  )
-
-               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
-               zckb    = ( cb0 + cb1 * zsqrt + cb2  * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr   &
-                  &    + ( cb5 + cb6 * zsqrt + cb7  * zsal )                                            &
-                  &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             &
-                  &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  )
-
-               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal
-               zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2
-
-               ! PKW (H2O) (DICKSON AND RILEY, 1979)
-               zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal
-
+               zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   &
+               &         + LOG(1.0d0 - 0.001005d0*zsal)            &
+               &         + LOG(1.0d0 + zst/zcks))
+
+
+               ! DISSOCIATION CONSTANT FOR BORATE
+               zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        &
+               &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         &
+               &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   &
+               &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      &
+               &      * zlogt + 0.053105*zsqrt*ztkel
+
+               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
+               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
+               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
+                  - 0.011555*zsal + 0.0001152*zsal*zsal)
+               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
+                  - 0.01781*zsal + 0.0001122*zsal*zsal)
+
+
+               ! PKW (H2O) (MILLERO, 1995) from composite data
+               zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    &
+                         - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
+
+               ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
+              zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   &
+              &          + (-106.736*ztr + 0.69171) * zsqrt       &
+              &          + (-0.65643*ztr - 0.01844) * zsal
+
+              zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
+              &          + (-160.340*ztr + 1.3566)*zsqrt          &
+              &          + (0.37335*ztr - 0.05778)*zsal
+
+              zck3p    = -3070.75*ztr - 18.126                    &
+              &          + (17.27039*ztr + 2.81197) * zsqrt       &
+              &          + (-44.99486*ztr - 0.09984) * zsal 
+
+              ! CONSTANT FOR SILICATE, MILLERO (1995)
+              zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
+              &          + (-458.79*ztr + 3.5913) * zisqrt       &
+              &          + (188.74*ztr - 1.5998) * zis           &
+              &          + (-12.1652*ztr + 0.07871) * zis2       &
+              &          + LOG(1.0 - 0.001005*zsal)
 
                ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
                !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
-               zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   &
-                  &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15
+               zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
+                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
+                  &      - 0.07711*zsal + 0.0041249*zsal15
+
+               ! CONVERT FROM DIFFERENT PH SCALES
+               total2free  = 1.0/(1.0 + zst/zcks)
+               free2SWS    = 1. + zst/zcks + zft/(zckf*total2free)
+               total2SWS   = total2free * free2SWS
+               SWS2total   = 1.0 / total2SWS
 
                ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
-               zak1    = 10**(zck1)
-               zak2    = 10**(zck2)
-               zakb    = EXP( zckb  )
+               zak1    = 10**(zck1) * total2SWS
+               zak2    = 10**(zck2) * total2SWS
+               zakb    = EXP( zckb ) * total2SWS
                zakw    = EXP( zckw )
                zaksp1  = 10**(zaksp0)
+               zak1p   = exp( zck1p )
+               zak2p   = exp( zck2p )
+               zak3p   = exp( zck3p )
+               zaksi   = exp( zcksi )
+               zckf    = zckf * total2SWS
 
                ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
@@ -327 +341 @@
                !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
                !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
-               zcpexp  = zpres /(rgas*ztkel)
-               zcpexp2 = zpres * zpres/(rgas*ztkel)
+               zcpexp  = zpres / (rgas*ztkel)
+               zcpexp2 = zpres * zcpexp
 
                ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
@@ -334 +348 @@
                !        (CF. BROECKER ET AL., 1982)
 
-               zbuf1  = -     ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
+               zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk40 + devk50 * ztc )
+               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
                zbuf2  = 0.5 * ( devk41 + devk51 * ztc )
-               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
                zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
                zbuf2  = 0.5 * ( devk42 + devk52 * ztc )
-               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
                zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
                zbuf2  = 0.5 * ( devk43 + devk53 * ztc )
-               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
                zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
                zbuf2  = 0.5 * ( devk44 + devk54 * ztc )
-               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
+               aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
+               akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk47 + devk57 * ztc )
+               ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk48 + devk58 * ztc )
+               ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk19 + devk29 * ztc + devk39 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk49 + devk59 * ztc )
+               ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk410 + devk510 * ztc )
+               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+               ! CONVERT FROM DIFFERENT PH SCALES
+               total2free  = 1.0/(1.0 + zst/aks3(ji,jj,jk))
+               free2SWS    = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk)
+               total2SWS   = total2free * free2SWS
+               SWS2total   = 1.0 / total2SWS
+
+               ! Convert to total scale
+               ak13(ji,jj,jk)  = ak13(ji,jj,jk)  * SWS2total
+               ak23(ji,jj,jk)  = ak23(ji,jj,jk)  * SWS2total
+               akb3(ji,jj,jk)  = akb3(ji,jj,jk)  * SWS2total
+               akw3(ji,jj,jk)  = akw3(ji,jj,jk)  * SWS2total
+               ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total
+               ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total
+               ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total
+               aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total
+               akf3(ji,jj,jk)  = akf3(ji,jj,jk)  / total2free
 
                ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
                !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
                !        (P. 1285) AND BERNER (1976)
-               zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
+               zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
+               zbuf2  = 0.5 * ( devk46 + devk56 * ztc )
                aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
 
-               ! TOTAL BORATE CONCENTR. [MOLES/L]
-               borat(ji,jj,jk) = bor1 * zcl * bor2
+               ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
+               borat(ji,jj,jk) = 0.0002414 * zcl / 10.811
+               sulfat(ji,jj,jk) = zst
+               fluorid(ji,jj,jk) = zft 
 
                ! Iron and SIO3 saturation concentration from ...
                sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
-               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) )
-
+               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel )
+
+               ! Liu and Millero (1999) only valid 5 - 50 degC
+               ztkel1 = MAX( 5. , tsn(ji,jj,jk,jp_tem) ) + 273.16
+               fesol(ji,jj,jk,1) = 10**((-13.486) - (0.1856* (zis**0.5)) + (0.3073*zis) + (5254/ztkel1))
+               fesol(ji,jj,jk,2) = 10**(2.517 - (0.885*(zis**0.5)) + (0.2139 * zis) - (1320/ztkel1) )
+               fesol(ji,jj,jk,3) = 10**(0.4511 - (0.3305*(ZIS**0.5)) - (1996/ztkel1) )
+               fesol(ji,jj,jk,4) = 10**(-0.2965 - (0.7881*(zis**0.5)) - (4086/ztkel1) )
+               fesol(ji,jj,jk,5) = 10**(4.4466 - (0.8505*(zis**0.5)) - (7980/ztkel1) )
             END DO
          END DO
@@ -373 +436 @@
    END SUBROUTINE p4z_che
 
+   SUBROUTINE ahini_for_at(p_hini)
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE ahini_for_at  ***
+      !!
+      !! Subroutine returns the root for the 2nd order approximation of the
+      !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic 
+      !! polynomial) around the local minimum, if it exists.
+      !! Returns * 1E-03_wp if p_alkcb <= 0
+      !!         * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot
+      !!         * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot
+      !!                    and the 2nd order approximation does not have 
+      !!                    a solution
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  ::  p_hini
+      INTEGER  ::   ji, jj, jk
+      REAL(wp)  ::  zca1, zba1
+      REAL(wp)  ::  zd, zsqrtd, zhmin
+      REAL(wp)  ::  za2, za1, za0
+      REAL(wp)  ::  p_dictot, p_bortot, p_alkcb 
+
+      IF( nn_timing == 1 )  CALL timing_start('ahini_for_at')
+      !
+      DO jk = 1, jpk
+        DO jj = 1, jpj
+          DO ji = 1, jpj
+            p_alkcb  = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn)
+            p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn)
+            p_bortot = borat(ji,jj,jk)
+            IF (p_alkcb <= 0.) THEN
+                p_hini(ji,jj,jk) = 1.e-3
+            ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN
+                p_hini(ji,jj,jk) = 1.e-10_wp
+            ELSE
+                zca1 = p_dictot/( p_alkcb + rtrn )
+                zba1 = p_bortot/ (p_alkcb + rtrn )
+           ! Coefficients of the cubic polynomial
+                za2 = aKb3(ji,jj,jk)*(1. - zba1) + ak13(ji,jj,jk)*(1.-zca1)
+                za1 = ak13(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - zca1)    &
+                &     + ak13(ji,jj,jk)*ak23(ji,jj,jk)*(1. - (zca1+zca1))
+                za0 = ak13(ji,jj,jk)*ak23(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1))
+                                        ! Taylor expansion around the minimum
+                zd = za2*za2 - 3.*za1   ! Discriminant of the quadratic equation
+                                        ! for the minimum close to the root
+
+                IF(zd > 0.) THEN        ! If the discriminant is positive
+                  zsqrtd = SQRT(zd)
+                  IF(za2 < 0) THEN
+                    zhmin = (-za2 + zsqrtd)/3.
+                  ELSE
+                    zhmin = -za1/(za2 + zsqrtd)
+                  ENDIF
+                  p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd)
+                ELSE
+                  p_hini(ji,jj,jk) = 1.e-7
+                ENDIF
+             !
+             ENDIF
+          END DO
+        END DO
+      END DO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('ahini_for_at')
+      !
+   END SUBROUTINE ahini_for_at
+
+   !===============================================================================
+   SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup )
+
+   ! Subroutine returns the lower and upper bounds of "non-water-selfionization"
+   ! contributions to total alkalinity (the infimum and the supremum), i.e
+   ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+])
+
+   ! Argument variables
+   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf
+   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup
+
+   p_alknw_inf(:,:,:) =  -trb(:,:,:,jppo4) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:)  &
+   &              - fluorid(:,:,:)
+   p_alknw_sup(:,:,:) =   (2. * trb(:,:,:,jpdic) + 2. * trb(:,:,:,jppo4) + trb(:,:,:,jpsil) )    &
+   &               * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:) 
+
+   END SUBROUTINE anw_infsup
+
+
+   SUBROUTINE solve_at_general( p_hini, zhi )
+
+   ! Universal pH solver that converges from any given initial value,
+   ! determines upper an lower bounds for the solution if required
+
+   ! Argument variables
+   !--------------------
+   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN)   :: p_hini
+   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  :: zhi
+
+   ! Local variables
+   !-----------------
+   INTEGER   ::  ji, jj, jk, jn
+   REAL(wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor
+   REAL(wp)  ::  zdelta, zh_delta
+   REAL(wp)  ::  zeqn, zdeqndh, zalka
+   REAL(wp)  ::  aphscale
+   REAL(wp)  ::  znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
+   REAL(wp)  ::  znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
+   REAL(wp)  ::  znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
+   REAL(wp)  ::  znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
+   REAL(wp)  ::  znumer_nh4, zdnumer_nh4, zdenom_nh4, zalk_nh4, zdalk_nh4
+   REAL(wp)  ::  znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s
+   REAL(wp)  ::  znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
+   REAL(wp)  ::  znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
+   REAL(wp)  ::  zalk_wat, zdalk_wat
+   REAL(wp)  ::  zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit
+   LOGICAL   ::  l_exitnow
+   REAL(wp), PARAMETER :: pz_exp_threshold = 1.0
+   REAL(wp), POINTER, DIMENSION(:,:,:) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
+
+   IF( nn_timing == 1 )  CALL timing_start('solve_at_general')
+      !  Allocate temporary workspace
+   CALL wrk_alloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask )
+   CALL wrk_alloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin )
+
+   CALL anw_infsup( zalknw_inf, zalknw_sup )
+
+   rmask(:,:,:) = tmask(:,:,:)
+   zhi(:,:,:)   = 0.
+
+   ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
+   DO jk = 1, jpk
+      DO jj = 1, jpj
+         DO ji = 1, jpj
+            IF (rmask(ji,jj,jk) == 1.) THEN
+               p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn)
+               aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
+               zh_ini = p_hini(ji,jj,jk)
+
+               zdelta = (p_alktot-zalknw_inf(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale
+
+               IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN
+                 zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) )
+               ELSE
+                 zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2.
+               ENDIF
+
+               zdelta = (p_alktot-zalknw_sup(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale
+
+               IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN
+                 zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2.
+               ELSE
+                 zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) )
+               ENDIF
+
+               zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk))
+            ENDIF
+         END DO
+      END DO
+   END DO
+
+   zeqn_absmin(:,:,:) = HUGE(1._wp)
+
+   DO jn = 1, jp_maxniter_atgen 
+   DO jk = 1, jpk
+      DO jj = 1, jpj
+         DO ji = 1, jpj
+            IF (rmask(ji,jj,jk) == 1.) THEN
+               zfact = rhop(ji,jj,jk) / 1000. + rtrn
+               p_alktot = trb(ji,jj,jk,jptal) / zfact
+               zdic  = trb(ji,jj,jk,jpdic) / zfact
+               zbot  = borat(ji,jj,jk)
+               zpt = trb(ji,jj,jk,jppo4) / zfact * po4r
+               zsit = trb(ji,jj,jk,jpsil) / zfact
+               zst = sulfat (ji,jj,jk)
+               zft = fluorid(ji,jj,jk)
+               aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
+               zh = zhi(ji,jj,jk)
+               zh_prev = zh
+
+               ! H2CO3 - HCO3 - CO3 : n=2, m=0
+               znumer_dic = 2.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)
+               zdenom_dic = ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*(ak13(ji,jj,jk) + zh)
+               zalk_dic   = zdic * (znumer_dic/zdenom_dic)
+               zdnumer_dic = ak13(ji,jj,jk)*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh     &
+                             *(4.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk))
+               zdalk_dic   = -zdic*(zdnumer_dic/zdenom_dic**2)
+
+
+               ! B(OH)3 - B(OH)4 : n=1, m=0
+               znumer_bor = akb3(ji,jj,jk)
+               zdenom_bor = akb3(ji,jj,jk) + zh
+               zalk_bor   = zbot * (znumer_bor/zdenom_bor)
+               zdnumer_bor = akb3(ji,jj,jk)
+               zdalk_bor   = -zbot*(zdnumer_bor/zdenom_bor**2)
+
+
+               ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
+               znumer_po4 = 3.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)  &
+               &            + zh*(2.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh* ak1p3(ji,jj,jk))
+               zdenom_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)     &
+               &            + zh*( ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh))
+               zalk_po4   = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1
+               zdnumer_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)  &
+               &             + zh*(4.*ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)         &
+               &             + zh*(9.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)                         &
+               &             + ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)                                &
+               &             + zh*(4.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh * ak1p3(ji,jj,jk) ) ) )
+               zdalk_po4   = -zpt * (zdnumer_po4/zdenom_po4**2)
+
+               ! H4SiO4 - H3SiO4 : n=1, m=0
+               znumer_sil = aksi3(ji,jj,jk)
+               zdenom_sil = aksi3(ji,jj,jk) + zh
+               zalk_sil   = zsit * (znumer_sil/zdenom_sil)
+               zdnumer_sil = aksi3(ji,jj,jk)
+               zdalk_sil   = -zsit * (zdnumer_sil/zdenom_sil**2)
+
+               ! HSO4 - SO4 : n=1, m=1
+               aphscale = 1.0 + zst/aks3(ji,jj,jk)
+               znumer_so4 = aks3(ji,jj,jk) * aphscale
+               zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh
+               zalk_so4   = zst * (znumer_so4/zdenom_so4 - 1.)
+               zdnumer_so4 = aks3(ji,jj,jk)
+               zdalk_so4   = -zst * (zdnumer_so4/zdenom_so4**2)
+
+               ! HF - F : n=1, m=1
+               znumer_flu =  akf3(ji,jj,jk)
+               zdenom_flu =  akf3(ji,jj,jk) + zh
+               zalk_flu   =  zft * (znumer_flu/zdenom_flu - 1.)
+               zdnumer_flu = akf3(ji,jj,jk)
+               zdalk_flu   = -zft * (zdnumer_flu/zdenom_flu**2)
+
+               ! H2O - OH
+               aphscale = 1.0 + zst/aks3(ji,jj,jk)
+               zalk_wat   = akw3(ji,jj,jk)/zh - zh/aphscale
+               zdalk_wat  = -akw3(ji,jj,jk)/zh**2 - 1./aphscale
+
+               ! CALCULATE [ALK]([CO3--], [HCO3-])
+               zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil   &
+               &      + zalk_so4 + zalk_flu                       &
+               &      + zalk_wat - p_alktot
+
+               zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil   &
+               &       + zalk_so4 + zalk_flu + zalk_wat)
+
+               zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil &
+               &         + zdalk_so4 + zdalk_flu + zdalk_wat
+
+               ! Adapt bracketing interval
+               IF(zeqn > 0._wp) THEN
+                 zh_min(ji,jj,jk) = zh_prev
+               ELSEIF(zeqn < 0._wp) THEN
+                 zh_max(ji,jj,jk) = zh_prev
+               ENDIF
+
+               IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN
+               ! if the function evaluation at the current point is
+               ! not decreasing faster than with a bisection step (at least linearly)
+               ! in absolute value take one bisection step on [ph_min, ph_max]
+               ! ph_new = (ph_min + ph_max)/2d0
+               !
+               ! In terms of [H]_new:
+               ! [H]_new = 10**(-ph_new)
+               !         = 10**(-(ph_min + ph_max)/2d0)
+               !         = SQRT(10**(-(ph_min + phmax)))
+               !         = SQRT(zh_max * zh_min)
+                  zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk))
+                  zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below
+               ELSE
+               ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH
+               !           = -zdeqndh * LOG(10) * [H]
+               ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10))
+               !
+               ! pH_new = pH_old + \deltapH
+               !
+               ! [H]_new = 10**(-pH_new)
+               !         = 10**(-pH_old - \Delta pH)
+               !         = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10)))
+               !         = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10)))
+               !         = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old))
+
+                  zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
+
+                  IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN
+                     zh          = zh_prev*EXP(zh_lnfactor)
+                  ELSE
+                     zh_delta    = zh_lnfactor*zh_prev
+                     zh          = zh_prev + zh_delta
+                  ENDIF
+
+                  IF( zh < zh_min(ji,jj,jk) ) THEN
+                     ! if [H]_new < [H]_min
+                     ! i.e., if ph_new > ph_max then
+                     ! take one bisection step on [ph_prev, ph_max]
+                     ! ph_new = (ph_prev + ph_max)/2d0
+                     ! In terms of [H]_new:
+                     ! [H]_new = 10**(-ph_new)
+                     !         = 10**(-(ph_prev + ph_max)/2d0)
+                     !         = SQRT(10**(-(ph_prev + phmax)))
+                     !         = SQRT([H]_old*10**(-ph_max))
+                     !         = SQRT([H]_old * zh_min)
+                     zh                = SQRT(zh_prev * zh_min(ji,jj,jk))
+                     zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below
+                  ENDIF
+
+                  IF( zh > zh_max(ji,jj,jk) ) THEN
+                     ! if [H]_new > [H]_max
+                     ! i.e., if ph_new < ph_min, then
+                     ! take one bisection step on [ph_min, ph_prev]
+                     ! ph_new = (ph_prev + ph_min)/2d0
+                     ! In terms of [H]_new:
+                     ! [H]_new = 10**(-ph_new)
+                     !         = 10**(-(ph_prev + ph_min)/2d0)
+                     !         = SQRT(10**(-(ph_prev + ph_min)))
+                     !         = SQRT([H]_old*10**(-ph_min))
+                     !         = SQRT([H]_old * zhmax)
+                     zh                = SQRT(zh_prev * zh_max(ji,jj,jk))
+                     zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below
+                  ENDIF
+               ENDIF
+
+               zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk))
+
+               ! Stop iterations once |\delta{[H]}/[H]| < rdel
+               ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel
+               ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
+
+               ! Alternatively:
+               ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))|
+               !             ~ 1/LOG(10) * |\Delta [H]|/[H]
+               !             < 1/LOG(10) * rdel
+
+               ! Hence |zeqn/(zdeqndh*zh)| < rdel
+
+               ! rdel <-- pp_rdel_ah_target
+               l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
+
+               IF(l_exitnow) THEN 
+                  rmask(ji,jj,jk) = 0.
+               ENDIF
+
+               zhi(ji,jj,jk) =  zh
+
+               IF(jn >= jp_maxniter_atgen) THEN
+                  zhi(ji,jj,jk) = -1._wp
+               ENDIF
+
+            ENDIF
+         END DO
+      END DO
+   END DO
+   END DO
+   !
+   CALL wrk_dealloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask )
+   CALL wrk_dealloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin )
+
+
+   IF( nn_timing == 1 )  CALL timing_stop('solve_at_general')
+
+
+   END SUBROUTINE solve_at_general
 
    INTEGER FUNCTION p4z_che_alloc()
@@ -378 +797 @@
       !!                     ***  ROUTINE p4z_che_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,2), chemo2(jpi,jpj,jpk), STAT=p4z_che_alloc )
+#if defined key_ligand
+      INTEGER ::   ierr(3)        ! Local variables
+#else
+      INTEGER ::   ierr(2)        ! Local variables
+#endif
+      !!----------------------------------------------------------------------
+
+      ierr(:) = 0
+
+      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) )
+
+      ALLOCATE( akb3(jpi,jpj,jpk)    ,                             &
+         &      akw3(jpi,jpj,jpk)    , borat (jpi,jpj,jpk) ,       &
+         &      aks3(jpi,jpj,jpk)    , akf3(jpi,jpj,jpk)   ,       &
+         &      ak1p3(jpi,jpj,jpk)   , ak2p3(jpi,jpj,jpk)  ,       &
+         &      ak3p3(jpi,jpj,jpk)   , aksi3(jpi,jpj,jpk)  ,       &
+         &      fluorid(jpi,jpj,jpk) , sulfat(jpi,jpj,jpk) ,     STAT=ierr(2) )
+
+      ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) )
+
+      !* Variable for chemistry of the CO2 cycle
+      p4z_che_alloc = MAXVAL( ierr )
       !
       IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.')
@@ -396 +836 @@
 
    !!======================================================================
-END MODULE p4zche
+END MODULE  p4zche

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

@@ -5 +5 @@
    !!======================================================================
    !! History :   3.5  !  2012-07 (O. Aumont, A. Tagliabue, C. Ethe) Original code
-   !!----------------------------------------------------------------------
-#if defined key_pisces
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
+   !!----------------------------------------------------------------------
+#if defined key_pisces || defined key_pisces_quota
    !!----------------------------------------------------------------------
    !!   'key_top'       and                                      TOP models
-   !!   'key_pisces'                                       PISCES bio-model
+   !!   'key_pisces*'                                       PISCES bio-model
    !!----------------------------------------------------------------------
    !!   p4z_fechem       :  Compute remineralization/scavenging of iron
@@ -33 +34 @@
    LOGICAL          ::  ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker
    LOGICAL          ::  ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
+   LOGICAL          ::  ln_fecolloid !: boolean for variable colloidal fraction
    REAL(wp), PUBLIC ::  xlam1        !: scavenging rate of Iron 
    REAL(wp), PUBLIC ::  xlamdust     !: scavenging rate of Iron by dust 
    REAL(wp), PUBLIC ::  ligand       !: ligand concentration in the ocean 
+   REAL(wp), PUBLIC ::  kfep         !: rate constant for nanoparticle formation
 
    REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth
@@ -48 +51 @@
 CONTAINS
 
-   SUBROUTINE p4z_fechem( kt, knt )
+   SUBROUTINE p4z_fechem( kt, jnt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_fechem  ***
@@ -62 +65 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
-      !
-      INTEGER  ::   ji, jj, jk, jic
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      !
+      INTEGER  ::   ji, jj, jk, jic, jn
       REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
-      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll
+      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, fe3sol
       REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag
       REAL(wp) ::   ztrc, zdust
@@ -72 +75 @@
       REAL(wp) ::   zdenom, zdenom2
 #endif
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig, precip
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zstrn, zstrn2
+      REAL(wp) ::   zzFeL1, zzFeL2, zzFe2, zzFeP, zzFe3, zzstrn2
+      REAL(wp) ::   zrum, zcodel, zargu, zval, zlight
       REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand
       REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2
       REAL(wp) :: zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2
-      REAL(wp) :: ztfe, zoxy
+      REAL(wp) :: ztfe, zoxy, zhplus
       REAL(wp) :: zstep
+#if defined key_ligand
+      REAL(wp) :: zaggliga, zaggligb
+      REAL(wp) :: dissol, zligco
+#endif
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -85 +95 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )
+      CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
       zFe3 (:,:,:) = 0.
       zFeL1(:,:,:) = 0.
       zTL1 (:,:,:) = 0.
       IF( ln_fechem ) THEN
+         CALL wrk_alloc( jpi, jpj,      zstrn, zstrn2 )
          CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
          zFe2 (:,:,:) = 0.
@@ -104 +115 @@
          ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
       ELSE
+#if defined key_ligand
+         ztotlig(:,:,:) = trb(:,:,:,jplgw) * 1E9
+#else
          ztotlig(:,:,:) = ligand * 1E9
+#endif
       ENDIF
 
       IF( ln_fechem ) THEN
+
+         ! compute the day length depending on latitude and the day
+         zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
+         zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
+
+         ! day length in hours
+         zstrn(:,:) = 0.
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
+               zargu = MAX( -1., MIN(  1., zargu ) )
+               zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
+            END DO
+         END DO
+
+         ! Maximum light intensity
+         zstrn2(:,:) = zstrn(:,:) / 24.
+         WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
+         zstrn(:,:) = 24. / zstrn(:,:)
+
          ! ------------------------------------------------------------
          ! NEW FE CHEMISTRY ROUTINE from Tagliabue and Volker (2009)
@@ -113 +148 @@
          ! Chemistry is supposed to be fast enough to be at equilibrium
          ! ------------------------------------------------------------
+         DO jn = 1, 2
 !CDIR NOVERRCHK
          DO jk = 1, jpkm1
@@ -119 +155 @@
 !CDIR NOVERRCHK
                DO ji = 1, jpi
+                  zlight = etot(ji,jj,jk) * zstrn(ji,jj) * float(2-jn)
+                  zzstrn2 = zstrn2(ji,jj) * float(2-jn) + (1. - zstrn2(ji,jj) ) * float(jn-1)
                   ! Calculate ligand concentrations : assume 2/3rd of excess goes to
                   ! strong ligands (L1) and 1/3rd to weak ligands (L2)
@@ -134 +172 @@
                   zkox   = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn )
                   ! PHOTOREDUCTION of complexed iron : Tagliabue and Arrigo (2006)
-                  zkph2 = MAX( 0., 15. * etot(ji,jj,jk) / ( etot(ji,jj,jk) + 2. ) )
+                  zkph2 = MAX( 0., 15. * zlight / ( zlight + 2. ) )
                   zkph1 = zkph2 / 5.
                   ! pass the dfe concentration from PISCES
@@ -180 +218 @@
                   ENDIF
                   ! solve for the other Fe species
-                  zFe3(ji,jj,jk) = MAX( 0., zxs ) 
-                  zFep(ji,jj,jk) = MAX( 0., ( ks * zFe3(ji,jj,jk) / kpr ) )
+                  zzFe3 = MAX( 0., zxs )
+                  zzFep = MAX( 0., ( ks * zzFe3 / kpr ) )
                   zkappa2 = ( kb2 + zkph2 ) / kl2
-                  zFeL2(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * zTL2(ji,jj,jk) ) / ( zkappa2 + zFe3(ji,jj,jk) ) )
-                  zFeL1(ji,jj,jk) = MAX( 0., ( ztfe / zb - za / zb * zFe3(ji,jj,jk) - zc / zb * zFeL2(ji,jj,jk) ) )
-                  zFe2 (ji,jj,jk) = MAX( 0., ( ( zkph1 * zFeL1(ji,jj,jk) + zkph2 * zFeL2(ji,jj,jk) ) / zkox ) )
+                  zzFeL2 = MAX( 0., ( zzFe3 * zTL2(ji,jj,jk) ) / ( zkappa2 + zzFe3 ) )
+                  zzFeL1 = MAX( 0., ( ztfe / zb - za / zb * zzFe3 - zc / zb * zzFeL2 ) )
+                  zzFe2  = MAX( 0., ( ( zkph1 * zzFeL1 + zkph2 * zzFeL2 ) / zkox ) )
+                  zFe3(ji,jj,jk)  = zFe3(ji,jj,jk)  + zzFe3 * zzstrn2
+                  zFe2(ji,jj,jk)  = zFe2(ji,jj,jk)  + zzFe2 * zzstrn2
+                  zFeL2(ji,jj,jk) = zFeL2(ji,jj,jk) + zzFeL2 * zzstrn2
+                  zFeL1(ji,jj,jk) = zFeL1(ji,jj,jk) + zzFeL1 * zzstrn2
+                  zFeP(ji,jj,jk)  = zFeP(ji,jj,jk)  + zzFeP * zzstrn2
                END DO
             END DO
+         END DO
          END DO
       ELSE
@@ -236 +280 @@
                   zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
                ELSE
-                  zfeequi = zFe3(ji,jj,jk) * 1E-9 
-                  zfecoll = 0.5 * zFeL1(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  &
+                     &         + 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 ) )
+                  ELSE
+                     zfeequi = zFe3(ji,jj,jk) * 1E-9 
+                     zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
+                     fe3sol  = 0.
+                     kfep    = 0.
+                  ENDIF
                ENDIF
 #if defined key_kriest
@@ -272 +327 @@
                    &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) + 5.09E3 * trb(ji,jj,jk,jppoc) )
                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
                !
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
+#  if defined key_ligand
+               zaggligb = zlam1b * zstep * zligco
+#  endif
+               ! precipitation of Fe3+, creation of nanoparticles
+               precip(ji,jj,jk) = max( 0., (zfeequi - fe3sol) ) * kfep * zstep
+               !
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb &
+               &                     - zcoag - precip(ji,jj,jk)
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
+#endif
+#if defined key_ligand
+               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
 #endif
             END DO
@@ -297 +374 @@
       ENDIF
 
+#if defined key_ligand
+      IF (.NOT. ln_fechem) THEN
+          plig(:,:,:) =  MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trn(:,:,:,jpfer) +rtrn ) ) )
+          plig(:,:,:) = max(0. , plig(:,:,:) )
+      ENDIF
+#endif
+
       !  Output of some diagnostics variables
       !     ---------------------------------
-      IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
-         IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
-         IF( iom_use("TL1")    )  CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
-         IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
-         IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! biron
-         IF( ln_fechem ) THEN
-            IF( iom_use("Fe2")  ) CALL iom_put("Fe2"    , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
-            IF( iom_use("FeL2") ) CALL iom_put("FeL2"   , zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
-            IF( iom_use("FeP")  ) CALL iom_put("FeP"    , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
-            IF( iom_use("TL2")  ) CALL iom_put("TL2"    , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
+      IF( ln_diatrc .AND. lk_iomput ) THEN
+         IF( jnt == nrdttrc ) THEN
+            CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
+            CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
+            CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
+            CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
+            CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! biron
+            IF( ln_fechem ) THEN
+               CALL iom_put("Fe2" , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
+               CALL iom_put("FeL2", zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
+               CALL iom_put("FeP" , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
+               CALL iom_put("TL2" , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
+            ENDIF
          ENDIF
       ENDIF
@@ -319 +405 @@
       ENDIF
       !
-                       CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )
+                       CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
       IF( ln_fechem )  CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
       !
@@ -339 +425 @@
       !!
       !!----------------------------------------------------------------------
-      NAMELIST/nampisfer/ ln_fechem, ln_ligvar, xlam1, xlamdust, ligand 
+      NAMELIST/nampisfer/ ln_fechem, ln_ligvar, ln_fecolloid, xlam1, xlamdust, ligand, kfep 
       INTEGER :: ios                 ! Local integer output status for namelist read
 
@@ -355 +441 @@
          WRITE(numout,*) ' Namelist parameters for Iron chemistry, nampisfer'
          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    enable complex iron chemistry scheme      ln_fechem =', ln_fechem
-         WRITE(numout,*) '    variable concentration of ligand          ln_ligvar =', ln_ligvar
-         WRITE(numout,*) '    scavenging rate of Iron                   xlam1     =', xlam1
-         WRITE(numout,*) '    scavenging rate of Iron by dust           xlamdust  =', xlamdust
-         WRITE(numout,*) '    ligand concentration in the ocean         ligand    =', ligand
+         WRITE(numout,*) '    enable complex iron chemistry scheme      ln_fechem    =', ln_fechem
+         WRITE(numout,*) '    variable concentration of ligand          ln_ligvar    =', ln_ligvar
+         WRITE(numout,*) '    Variable colloidal fraction of Fe3+       ln_fecolloid =', ln_fecolloid
+         WRITE(numout,*) '    scavenging rate of Iron                   xlam1        =', xlam1
+         WRITE(numout,*) '    scavenging rate of Iron by dust           xlamdust     =', xlamdust
+         WRITE(numout,*) '    ligand concentration in the ocean         ligand       =', ligand
+         WRITE(numout,*) '    rate constant for nanoparticle formation  kfep         =', kfep
       ENDIF
       !

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

@@ -10 +10 @@
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
    !!                  !  2011-02  (J. Simeon, J. Orr) Include total atm P correction 
-   !!----------------------------------------------------------------------
-#if defined key_pisces
-   !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
+   !!----------------------------------------------------------------------
+#if defined key_pisces || defined key_pisces_quota
+   !!----------------------------------------------------------------------
+   !!   'key_pisces*'                                      PISCES bio-model
    !!----------------------------------------------------------------------
    !!   p4z_flx       :   CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
@@ -52 +53 @@
    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:)  ::  patm      ! atmospheric pressure at kt                 [N/m2]
    TYPE(FLD), ALLOCATABLE,       DIMENSION(:)    ::  sf_patm   ! structure of input fields (file informations, fields read)
-
+   LOGICAL, PUBLIC ::   ln_presatmco2  !: accounting for spatial atm CO2 in the compuation of carbon flux (T) or not (F)
+   TYPE(FLD), ALLOCATABLE,       DIMENSION(:)    ::  sf_atmco2   ! structure of input fields (file informations, fields read)
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: oce_co2   !: ocean carbon flux 
@@ -84 +86 @@
       !
       INTEGER  ::   ji, jj, jm, iind, iindm1
-      REAL(wp) ::   ztc, ztc2, ztc3, zws, zkgwan
+      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan
       REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact
-      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2
+      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff
+      REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2
       REAL(wp) ::   zyr_dec, zdco2dt
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d 
+      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d
       !!---------------------------------------------------------------------
       !
@@ -103 +106 @@
       IF( kt /= nit000 .AND. knt == 1 ) CALL p4z_patm( kt )    ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
 
-      IF( ln_co2int ) THEN 
+      IF( ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
          ! Linear temporal interpolation  of atmospheric pco2.  atcco2.txt has annual values.
          ! Caveats: First column of .txt must be in years, decimal  years preferably. 
@@ -121 +124 @@
 #endif
 
-      DO jm = 1, 10
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-
-               ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
-               zbot  = borat(ji,jj,1)
-               zfact = rhop(ji,jj,1) / 1000. + rtrn
-               zdic  = trb(ji,jj,1,jpdic) / zfact
-               zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
-               zalka = trb(ji,jj,1,jptal) / zfact
-
-               ! CALCULATE [ALK]([CO3--], [HCO3-])
-               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph + zbot / ( 1.+ zph / akb3(ji,jj,1) )  )
-
-               ! CALCULATE [H+] AND [H2CO3]
-               zah2   = SQRT(  (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1)   &
-                  &                                        / ak13(ji,jj,1) ) * ( 2.* zdic - zalk )  )
-               zah2   = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 )
-               zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact
-               hi(ji,jj,1)   = zah2 * zfact
-            END DO
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
+            zfact = rhop(ji,jj,1) / 1000. + rtrn
+            zdic  = trb(ji,jj,1,jpdic) / zfact
+            zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
+            ! CALCULATE [H2CO3]
+            zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)*zfact
          END DO
       END DO
-
 
       ! --------------
@@ -162 +149 @@
             ztc2 = ztc * ztc
             ztc3 = ztc * ztc2 
+            ztc4 = ztc2 * ztc2
             ! Compute the schmidt Number both O2 and CO2
-            zsch_co2 = 2073.1 - 125.62 * ztc + 3.6276 * ztc2 - 0.043126 * ztc3
-            zsch_o2  = 1953.4 - 128.0  * ztc + 3.9918 * ztc2 - 0.050091 * ztc3
+            zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4
+            zsch_o2  = 1920.4 - 135.6  * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4
             !  wind speed 
             zws  = wndm(ji,jj) * wndm(ji,jj)
             ! Compute the piston velocity for O2 and CO2
-            zkgwan = 0.3 * zws  + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946  * ztc2 )
+            zkgwan = 0.251 * zws
             zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
 # if defined key_degrad
@@ -181 +169 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
+            ztkel = tsn(ji,jj,1,jp_tem) + 273.15
+            zsal  = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35.
+            zvapsw = exp(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*log(ztkel/100) - 0.000544*zsal)
+            xCO2approx = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) * 1.E-6
+            zxc2 = (1.0 - xCO2approx)**2
+            zfugcoeff = exp(patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   &
+            &           / (82.05736 * ztkel))
+            zfco2 = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) * zfugcoeff
             ! Compute CO2 flux for the sea and air
-            zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)   ! (mol/L) * (m/s)
+            zfld = zfco2 * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)   ! (mol/L) * (m/s)
             zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
             oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000.
             ! compute the trend
             tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1)
-
             ! Compute O2 flux 
-            zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
+            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
             zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
             zoflx(ji,jj) = zfld16 - zflu16
@@ -198 +193 @@
       t_oce_co2_flx     = glob_sum( oce_co2(:,:) )                    !  Total Flux of Carbon
       t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx       !  Cumulative Total Flux of Carbon
-!      t_atm_co2_flx     = glob_sum( satmco2(:,:) * e1e2t(:,:) )       ! Total atmospheric pCO2
-      t_atm_co2_flx     =  atcco2      ! Total atmospheric pCO2
- 
+      t_atm_co2_flx = glob_sum( satmco2(:,:) * e1e2t(:,:) )         ! Total atmospheric pCO2
+
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('flx ')")
@@ -208 +202 @@
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         CALL wrk_alloc( jpi, jpj, zw2d )  
+         CALL wrk_alloc( jpi, jpj, zw2d )
          IF( iom_use( "Cflx"  ) )  THEN
             zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r
-            CALL iom_put( "Cflx"     , zw2d ) 
+            CALL iom_put( "Cflx"     , zw2d )
          ENDIF
          IF( iom_use( "Oflx"  ) )  THEN
@@ -222 +216 @@
          ENDIF
          IF( iom_use( "Dpco2" ) ) THEN
-           zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
+           zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) / 1000. * rfact2r * tmask(:,:,1) / ( zkgco2(:,:) * chemc(:,:,1) + rtrn )
            CALL iom_put( "Dpco2" ,  zw2d )
          ENDIF
          IF( iom_use( "Dpo2" ) )  THEN
-           zw2d(:,:) = ( atcox * patm(:,:) - trb(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1)
+           zw2d(:,:) = ( patm(:,:) - trb(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * atcox * tmask(:,:,1)
            CALL iom_put( "Dpo2"  , zw2d )
          ENDIF
@@ -235 +229 @@
       ELSE
          IF( ln_diatrc ) THEN
-            trc2d(:,:,jp_pcs0_2d    ) = oce_co2(:,:) / e1e2t(:,:) * rfact2r 
+            trc2d(:,:,jp_pcs0_2d    ) = oce_co2(:,:) / e1e2t(:,:) * rfact2r
+            trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1)
+            trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1)
+            trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
+         ENDIF
+      ENDIF
+
+      IF( ln_diatrc ) THEN
+         IF( lk_iomput .AND. knt == nrdttrc ) THEN
+            CALL iom_put( "Cflx" , oce_co2(:,:) / e1e2t(:,:) / rfact ) 
+            CALL iom_put( "Oflx" , zoflx(:,:) * 1000 * tmask(:,:,1)  )
+            CALL iom_put( "Kg"   , zkgco2(:,:) * tmask(:,:,1) )
+            CALL iom_put( "Dpco2", ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) )
+            CALL iom_put( "Dpo2" , ( atcox * patm(:,:) - trb(:,:,1,jpoxy) * atcox / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
+         ELSE
+            trc2d(:,:,jp_pcs0_2d    ) = oce_co2(:,:) / e1e2t(:,:) / rfact 
             trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 
             trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) 
@@ -281 +290 @@
          WRITE(numout,*) ' '
       ENDIF
-      IF( .NOT.ln_co2int ) THEN
+      IF( .NOT.ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
          IF(lwp) THEN                         ! control print
             WRITE(numout,*) '    Constant Atmospheric pCO2 value  atcco2    =', atcco2
@@ -287 +296 @@
          ENDIF
          satmco2(:,:)  = atcco2      ! Initialisation of atmospheric pco2
-      ELSE
+      ELSEIF(ln_co2int .AND. .NOT.ln_presatmco2) THEN
          IF(lwp)  THEN
             WRITE(numout,*) '    Atmospheric pCO2 value  from file clname      =', TRIM( clname )
@@ -309 +318 @@
          END DO
          CLOSE(numco2)
+      ELSEIF(.NOT.ln_co2int .AND. ln_presatmco2) THEN
+         IF(lwp)  THEN
+            WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
+            WRITE(numout,*) ' '
+         ENDIF
+      ELSE
+         IF(lwp)  THEN
+            WRITE(numout,*) '    Spatialized Atmospheric pCO2 from an external file'
+            WRITE(numout,*) ' '
+         ENDIF
       ENDIF
       !
       oce_co2(:,:)  = 0._wp                ! Initialization of Flux of Carbon
+      t_atm_co2_flx = 0._wp
       t_oce_co2_flx = 0._wp
-      t_atm_co2_flx = 0._wp
       !
       CALL p4z_patm( nit000 )
@@ -335 +354 @@
       CHARACTER(len=100) ::  cn_dir   ! Root directory for location of ssr files
       TYPE(FLD_N)        ::  sn_patm  ! informations about the fields to be read
-      !!
-      NAMELIST/nampisatm/ ln_presatm, sn_patm, cn_dir
-
+      TYPE(FLD_N)        ::  sn_atmco2 ! informations about the fields to be read
+      !!
+      NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir
       !                                         ! ----------------------- !
       IF( kt == nit000 ) THEN                   ! First call kt=nittrc000 !
@@ -355 +374 @@
             WRITE(numout,*) '   Namelist nampisatm : Atmospheric Pressure as external forcing'
             WRITE(numout,*) '      constant atmopsheric pressure (F) or from a file (T)  ln_presatm = ', ln_presatm
+            WRITE(numout,*) '      spatial atmopsheric CO2 for flux calcs  ln_presatmco2 = ', ln_presatmco2
             WRITE(numout,*)
          ENDIF
@@ -367 +387 @@
          ENDIF
          !                                         
-         IF( .NOT.ln_presatm )   patm(:,:) = 1.e0    ! Initialize patm if no reading from a file
+         IF( ln_presatmco2 ) THEN
+            ALLOCATE( sf_atmco2(1), STAT=ierr )           !* allocate and fill sf_atmco2 (forcing structure) with sn_atmco2
+            IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_atmco2 structure' )
+            !
+            CALL fld_fill( sf_atmco2, (/ sn_atmco2 /), cn_dir, 'p4z_flx', 'Atmospheric co2 partial pressure ', 'nampisatm' )
+                                   ALLOCATE( sf_atmco2(1)%fnow(jpi,jpj,1)   )
+            IF( sn_atmco2%ln_tint )  ALLOCATE( sf_atmco2(1)%fdta(jpi,jpj,1,2) )
+         ENDIF
          !
       ENDIF
@@ -374 +401 @@
          CALL fld_read( kt, 1, sf_patm )               !* input Patm provided at kt + 1/2
          patm(:,:) = sf_patm(1)%fnow(:,:,1)                        ! atmospheric pressure
+      ELSE
+         patm(:,:) = 1.e0    ! Initialize patm if no reading from a file
+      ENDIF
+      !
+      IF( ln_presatmco2 ) THEN
+         CALL fld_read( kt, 1, sf_atmco2 )               !* input atmco2 provided at kt + 1/2
+         satmco2(:,:) = sf_atmco2(1)%fnow(:,:,1)                        ! atmospheric pressure
+      ELSE
+         satmco2(:,:) = atcco2    ! Initialize atmco2 if no reading from a file
+      ENDIF
+      !
+      IF(lwp) THEN                                 !* control print
+         WRITE(numout,*) ' Min-Max atmospheric CO2 = ', MINVAL(satmco2(:,:)),MAXVAL(satmco2(:,:))
       ENDIF
       !
@@ -400 +440 @@
 
    !!======================================================================
-END MODULE p4zflx
+END MODULE  p4zflx

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

@@ -6 +6 @@
    !! History :   1.0  !  2004-03 (O. Aumont) Original code
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
-#if defined key_pisces
+#if defined key_pisces || defined key_pisces_quota
    !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
+   !!   'key_pisces*'                                       PISCES bio-model
    !!----------------------------------------------------------------------
    !!   p4z_int        :  interpolation and computation of various accessory fields
@@ -81 +82 @@
 
    !!======================================================================
-END MODULE p4zint
+END MODULE  p4zint

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

@@ -26 +26 @@
    PUBLIC p4z_lim    
    PUBLIC p4z_lim_init    
+   PUBLIC p4z_lim_alloc
 
    !! * Shared module variables
@@ -47 +48 @@
    REAL(wp), PUBLIC ::  qdfelim     !:  optimal Fe quota for diatoms
    REAL(wp), PUBLIC ::  caco3r      !:  mean rainratio 
+
+   !!* Phytoplankton limitation terms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanono3   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatno3   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanonh4   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatnh4   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanopo4   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatpo4   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimphy    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdia    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimnfe    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdfe    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimsi     !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimbac    !: ??
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimbacl   !: ??
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concdfe    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concnfe    !: ???
 
    ! Coefficient for iron limitation
@@ -255 +273 @@
    END SUBROUTINE p4z_lim_init
 
+   INTEGER FUNCTION p4z_lim_alloc()
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE p5z_lim_alloc  ***
+      !!----------------------------------------------------------------------
+      USE lib_mpp , ONLY: ctl_warn
+      !!----------------------------------------------------------------------
+
+      !*  Biological arrays for phytoplankton growth
+      ALLOCATE( xnanono3(jpi,jpj,jpk), xdiatno3(jpi,jpj,jpk),       &
+         &      xnanonh4(jpi,jpj,jpk), xdiatnh4(jpi,jpj,jpk),       &
+         &      xnanopo4(jpi,jpj,jpk), xdiatpo4(jpi,jpj,jpk),       &
+         &      xlimphy (jpi,jpj,jpk), xlimdia (jpi,jpj,jpk),       &
+         &      xlimnfe (jpi,jpj,jpk), xlimdfe (jpi,jpj,jpk),       &
+         &      xlimbac (jpi,jpj,jpk), xlimbacl(jpi,jpj,jpk),       &
+         &      concnfe (jpi,jpj,jpk), concdfe (jpi,jpj,jpk),       &
+         &      xlimsi  (jpi,jpj,jpk), STAT=p4z_lim_alloc )
+      !
+      IF( p4z_lim_alloc /= 0 ) CALL ctl_warn('p4z_lim_alloc : failed to allocate arrays.')
+      !
+   END FUNCTION p4z_lim_alloc
+
+
 #else
    !!======================================================================

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

@@ -11 +11 @@
    !!                  !  2011-02  (J. Simeon, J. Orr)  Calcon salinity dependence
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Improvment of calcite dissolution
-   !!----------------------------------------------------------------------
-#if defined key_pisces
-   !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
+   !!----------------------------------------------------------------------
+#if defined key_pisces || defined key_pisces_quota
+   !!----------------------------------------------------------------------
+   !!   'key_pisces*'                                       PISCES bio-model
    !!----------------------------------------------------------------------
    !!   p4z_lys        :   Compute the CaCO3 dissolution 
@@ -22 +23 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
+   USE p4zche          !  Chemical model
    USE prtctl_trc      !  print control for debugging
    USE iom             !  I/O manager
@@ -60 +62 @@
       !
       INTEGER, INTENT(in) ::   kt, knt ! ocean time step
-      INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zalk, zdic, zph, zah2
-      REAL(wp) ::   zdispot, zfact, zcalcon, zalka, zaldi
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zdispot, zfact, zcalcon
       REAL(wp) ::   zomegaca, zexcess, zexcess0
+      REAL(wp) ::   zrfact2
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss   
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss, zhinit, zhi
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_lys')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss )
+      CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi )
       !
       zco3    (:,:,:) = 0.
       zcaldiss(:,:,:) = 0.
+      zhinit(:,:,:)   = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
       !     -------------------------------------------
       !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
       !     -------------------------------------------
       
-      DO jn = 1, 5                               !  BEGIN OF ITERATION
-         !
-!CDIR NOVERRCHK
-         DO jk = 1, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  zfact = rhop(ji,jj,jk) / 1000. + rtrn
-                  zph  = hi(ji,jj,jk) * tmask(ji,jj,jk) / zfact + ( 1.-tmask(ji,jj,jk) ) * 1.e-9 ! [H+]
-                  zdic  = trb(ji,jj,jk,jpdic) / zfact
-                  zalka = trb(ji,jj,jk,jptal) / zfact
-                  ! CALCULATE [ALK]([CO3--], [HCO3-])
-                  zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) )
-                  ! CALCULATE [H+] and [CO3--]
-                  zaldi = zdic - zalk
-                  zah2  = SQRT( zaldi * zaldi + 4.* ( zalk * ak23(ji,jj,jk) / ak13(ji,jj,jk) ) * ( zdic + zaldi ) )
-                  zah2  = 0.5 * ak13(ji,jj,jk) / zalk * ( zaldi + zah2 )
-                  !
-                  zco3(ji,jj,jk) = zalk / ( 2. + zah2 / ak23(ji,jj,jk) ) * zfact
-                  hi(ji,jj,jk)   = zah2 * zfact
-               END DO
+      CALL solve_at_general(zhinit, zhi)
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zfact = rhop(ji,jj,jk) / 1000. + rtrn
+               zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
+               &                + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
+               hi(ji,jj,jk)   = zhi(ji,jj,jk) * zfact
             END DO
          END DO
-         !
       END DO 
 
@@ -125 +114 @@
                zexcess0 = MAX( 0., excess(ji,jj,jk) )
                zexcess  = zexcess0**nca
-
                ! AMOUNT CACO3 (12C) THAT RE-ENTERS SOLUTION
                !       (ACCORDING TO THIS FORMULATION ALSO SOME PARTICULATE
@@ -133 +121 @@
                zdispot = zdispot * facvol(ji,jj,jk)
 # endif
-              !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
-              !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
-              zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
-              zco3(ji,jj,jk)      = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk)
-              !
-              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk)
-              tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) -      zcaldiss(ji,jj,jk)
-              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +      zcaldiss(ji,jj,jk)
+               !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
+               !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
+               zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss  ! calcite dissolution
+               !
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk)
+               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) -      zcaldiss(ji,jj,jk)
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +      zcaldiss(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-
+         !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
          IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )          * tmask(:,:,:) )
@@ -151 +138 @@
          IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", aksp(:,:,:) * 1.e+3 / calcon      * tmask(:,:,:) )
          IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) )
+
       ELSE
          IF( ln_diatrc ) THEN
@@ -158 +146 @@
          ENDIF
       ENDIF
+         !
       !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
@@ -165 +154 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss )
+      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_lys')
@@ -184 +173 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER  ::  ji, jj, jk
       INTEGER  ::  ios                 ! Local integer output status for namelist read
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
 
       NAMELIST/nampiscal/ kdca, nca
@@ -225 +211 @@
 #endif 
    !!======================================================================
-END MODULE p4zlys
+END MODULE  p4zlys

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

@@ -73 +73 @@
       REAL(wp) :: zgraze2 , zdenom, zdenom2
       REAL(wp) :: zfact   , zstep, zfood, zfoodlim, zproport
-      REAL(wp) :: zmortzgoc, zfrac, zfracfe, zratio, zratio2
+      REAL(wp) :: zmortzgoc, zfrac, zfracfe, zratio, zratio2, zfracal, zgrazcal
       REAL(wp) :: zepshert, zepsherv, zgrarsig, zgraztot, zgraztotn, zgraztotf
       REAL(wp) :: zgrarem2, zgrafer2, zgrapoc2, zprcaca, zmortz2, zgrasrat, zgrasratn
@@ -204 +204 @@
                tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgrarsig
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem2 - zgrarsig
+#if defined key_ligand
+               tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem2 - zgrarsig) * ldocz
+#endif
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarsig
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer2
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarsig
                tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zgrarsig              
-
                zmortz2 = ztortz2 + zrespz2
                zmortzgoc = unass2 / ( 1. - epsher2 ) * ztortz2 + zrespz2
@@ -222 +224 @@
                tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazf
 
-               ! calcite production
-               zprcaca = xfracal(ji,jj,jk) * zgrazn
-               prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
-               !
-               zprcaca = part2 * zprcaca
-               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * zprcaca
-               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca
 #if defined key_kriest
               znumpoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )
@@ -237 +231 @@
               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortzgoc - zgrazfffp - zgrazpof    &
                  &                 + zgraztotf * unass2
+              zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + rtrn )
+              zgrazcal = ( zgrazffep + zgrazpoc ) * (1. - part2) * zfracal
 #else
               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zgrazpoc - zgrazffep + zfrac
+              prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zfrac
+              conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazpoc - zgrazffep
               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zmortzgoc - zgrazffeg + zgrapoc2 - zfrac
+              prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zmortzgoc + zgrapoc2
+              consgoc(ji,jj,jk) = consgoc(ji,jj,jk) - zgrazffeg - zfrac
               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zgrazpof - zgrazfffp + zfracfe
               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + ferat3 * zmortzgoc - zgrazfffg     &
                  &                + zgraztotf * unass2 - zfracfe
+              zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + rtrn )
+              zgrazcal = zgrazffeg * (1. - part2) * zfracal
 #endif
+
+               ! calcite production
+               zprcaca = xfracal(ji,jj,jk) * zgrazn
+               prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
+               !
+               zprcaca = part2 * zprcaca
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrazcal - zprcaca
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * ( zgrazcal + zprcaca )
+               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) - zgrazcal + zprcaca
+
             END DO
          END DO

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

@@ -20 +20 @@
    USE p4zlim          !  Co-limitations
    USE p4zsink         !  vertical flux of particulate matter due to sinking
-   USE p4zint          !  interpolation and computation of various fields
    USE p4zprod         !  production
    USE iom             !  I/O manager
@@ -150 +149 @@
                tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgrarsig
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem - zgrarsig
+#if defined key_ligand
+               tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem - zgrarsig) * ldocz
+#endif
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarsig
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer
                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zgrapoc
+               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zgrapoc
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zgraztotf * unass
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarsig
@@ -172 +175 @@
                tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazsf
                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz - zgrazm
+               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zmortz
+               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazm
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz - zgrazmf
                !

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

@@ -17 +17 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p4zsink         !  vertical flux of particulate matter due to sinking
    USE p4zprod         !  Primary productivity 
+   USE p4zlim          !  Phytoplankton limitation terms
    USE prtctl_trc      !  print control for debugging
 
@@ -128 +128 @@
                tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfracal * zmortp
                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + ( 1. - zfracal ) * zmortp
+               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + ( 1. - zfracal ) * zmortp
+               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zfracal * zmortp
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ( 1. - zfracal ) * zmortp * zfactfe
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zfracal * zmortp * zfactfe
@@ -211 +213 @@
                tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zrespp2 + 0.5 * ztortp2
                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + 0.5 * ztortp2
+               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + 0.5 * ztortp2
+               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zrespp2 + 0.5 * ztortp2
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 0.5 * ztortp2 * zfactfe
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + ( zrespp2 + 0.5 * ztortp2 ) * zfactfe

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

@@ -9 +9 @@
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Improve light availability of nano & diat
    !!----------------------------------------------------------------------
-#if defined  key_pisces
-   !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
+#if defined  key_pisces || defined key_pisces_quota
+   !!----------------------------------------------------------------------
+   !!   'key_pisces*'                                      PISCES bio-model
    !!----------------------------------------------------------------------
    !!   p4z_opt       : light availability in the water column
@@ -43 +43 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: enano, ediat   !: PAR for phyto, nano and diat 
+#if defined key_pisces_quota
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: epico          !: PAR for pico
+#endif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etot_ndcy      !: PAR over 24h in case of diurnal cycle
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emoy           !: averaged PAR in the mixed layer
@@ -77 +80 @@
       REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100
+#if defined key_pisces_quota
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zetmp5
+#endif
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3
       !!---------------------------------------------------------------------
@@ -85 +91 @@
       CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
       CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 )
+#if defined key_pisces_quota
+      CALL wrk_alloc( jpi, jpj,      zetmp5 )
+#endif
 
       IF( knt == 1 .AND. ln_varpar ) CALL p4z_opt_sbc( kt )
@@ -99 +108 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
+#if defined key_pisces_quota
+               zchl = ( trb(ji,jj,jk,jpnch) + trb(ji,jj,jk,jpdch) + trb(ji,jj,jk,jppch) + rtrn ) * 1.e6
+#else
                zchl = ( trb(ji,jj,jk,jpnch) + trb(ji,jj,jk,jpdch) + rtrn ) * 1.e6
+#endif
                zchl = MIN(  10. , MAX( 0.05, zchl )  )
                irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
@@ -121 +134 @@
             enano    (:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
             ediat    (:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk)
+#if defined key_pisces_quota
+            epico    (:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
+#endif
          END DO
          !
@@ -139 +155 @@
             enano(:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
             ediat(:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk)
+#if defined key_pisces_quota
+            epico(:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
+#endif
          END DO
          etot_ndcy(:,:,:) =  etot(:,:,:) 
@@ -177 +196 @@
       zetmp3 (:,:)   = 0.e0
       zetmp4 (:,:)   = 0.e0
+#if defined key_pisces_quota
+      zetmp5 (:,:)   = 0.e0
+#endif
 
       DO jk = 1, nksrp
@@ -188 +210 @@
                   zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
                   zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
+#if defined key_pisces_quota
+                  zetmp5 (ji,jj) = zetmp5 (ji,jj) + epico    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
+#endif
                   zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk)
                ENDIF
@@ -208 +233 @@
                   enano(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
                   ediat(ji,jj,jk) = zetmp4(ji,jj) * z1_dep
+#if defined key_pisces_quota
+                  epico(ji,jj,jk) = zetmp5(ji,jj) * z1_dep
+#endif
                ENDIF
             END DO
@@ -228 +256 @@
       CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
       CALL wrk_dealloc( jpi, jpj, jpk, zpar,  ze0, ze1, ze2, ze3 )
+#if defined key_pisces_quota
+      CALL wrk_dealloc( jpi, jpj,      zetmp5 )
+#endif
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_opt')
@@ -410 +441 @@
                          enano    (:,:,:) = 0._wp
                          ediat    (:,:,:) = 0._wp
+#if defined key_pisces_quota
+                         epico    (:,:,:) = 0._wp
+#endif
       IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp
       ! 
@@ -423 +457 @@
       ALLOCATE( ekb(jpi,jpj,jpk)      , ekr(jpi,jpj,jpk), ekg(jpi,jpj,jpk),   &
         &       enano(jpi,jpj,jpk)    , ediat(jpi,jpj,jpk), &
+#if defined key_pisces_quota
+        &       epico(jpi,jpj,jpk)    ,                     &
+#endif
         &       etot_ndcy(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc ) 
          !

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

@@ -80 +80 @@
       REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap
       REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zproreg, zproreg2
-      REAL(wp) ::   zmxltst, zmxlday, zmaxday
+      REAL(wp) ::   zmxltst, zmxlday, zmaxday, zdocprod
       REAL(wp) ::   zpislopen  , zpislope2n
-      REAL(wp) ::   zrum, zcodel, zargu, zval
+      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
       REAL(wp) ::   zfact
       CHARACTER (len=25) :: charout
@@ -390 +390 @@
               zproreg  = zprorca(ji,jj,jk) - zpronew(ji,jj,jk)
               zproreg2 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk)
+              zdocprod = excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk)
               tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronew(ji,jj,jk) - zpronewd(ji,jj,jk)
@@ -400 +401 @@
               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcret2
               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk)
+              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zdocprod
+              zfeup = texcret * zprofen(ji,jj,jk) + texcret2 * zprofed(ji,jj,jk)
+#if defined key_ligand
+              tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp  - zfeup * plig(ji,jj,jk) * lthet
+#endif
               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
                  &                + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
-              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - texcret * zprofen(ji,jj,jk) - texcret2 * zprofed(ji,jj,jk)
+              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup
               tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)

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

@@ -7 +7 @@
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Quota model for iron
+   !!             3.6  !  2016-03  (O. Aumont) Quota model and reorganization
    !!----------------------------------------------------------------------
 #if defined key_pisces
@@ -22 +23 @@
    USE p4zopt          !  optical model
    USE p4zche          !  chemical model
+   USE p4zlim          ! Phytoplankton limitation factors
    USE p4zprod         !  Growth rate of the 2 phyto groups
-   USE p4zmeso         !  Sources and sinks of mesozooplankton
-   USE p4zint          !  interpolation and computation of various fields
-   USE p4zlim
    USE prtctl_trc      !  print control for debugging
    USE iom             !  I/O manager
@@ -38 +37 @@
 
    !! * Shared module variables
-   REAL(wp), PUBLIC ::  xremik     !: remineralisation rate of POC 
-   REAL(wp), PUBLIC ::  xremip     !: remineralisation rate of DOC
+   REAL(wp), PUBLIC ::  xremik     !: remineralisation rate of DOC
    REAL(wp), PUBLIC ::  nitrif     !: NH4 nitrification rate 
-   REAL(wp), PUBLIC ::  xsirem     !: remineralisation rate of POC 
-   REAL(wp), PUBLIC ::  xsiremlab  !: fast remineralisation rate of POC 
+   REAL(wp), PUBLIC ::  xsirem     !: remineralisation rate of BSi
+   REAL(wp), PUBLIC ::  xsiremlab  !: fast remineralisation rate of BSi
    REAL(wp), PUBLIC ::  xsilab     !: fraction of labile biogenic silica 
    REAL(wp), PUBLIC ::  oxymin     !: halk saturation constant for anoxia 
 
-
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitr     !: denitrification array
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitnh4   !: -    -    -    -   -
+
 
    !!* Substitution
@@ -59 +57 @@
 CONTAINS
 
-   SUBROUTINE p4z_rem( kt, knt )
+   SUBROUTINE p4z_rem( kt, jnt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_rem  ***
@@ -68 +66 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zremip, zremik, zsiremin 
+      REAL(wp) ::   zremik, zsiremin 
       REAL(wp) ::   zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
-      REAL(wp) ::   zbactfer, zorem, zorem2, zofer, zolimit
+      REAL(wp) ::   zbactfer, zolimit
       REAL(wp) ::   zosil, ztem
-#if ! defined key_kriest
-      REAL(wp) ::   zofer2
-#endif
-      REAL(wp) ::   zonitr, zstep, zfact
+      REAL(wp) ::   zonitr, zstep, zrfact2
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztempbac
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod, zw3d
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod
       !!---------------------------------------------------------------------
       !
@@ -87 +82 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      ztempbac                  )
+      CALL wrk_alloc( jpi, jpj,      ztempbac )
       CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi )
+
+      ! Initialization of local variables
+      ! ---------------------------------
 
       ! Initialisation of temprary arrys
@@ -219 +217 @@
                zstep = zstep * facvol(ji,jj,jk)
 # endif
-               ! POC disaggregation by turbulence and bacterial activity. 
-               ! --------------------------------------------------------
-               zremip = xremip * zstep * tgfunc(ji,jj,jk) * ( 1.- 0.55 * nitrfac(ji,jj,jk) ) 
-
-               ! POC disaggregation rate is reduced in anoxic zone as shown by
-               ! sediment traps data. In oxic area, the exponent of the martin s
-               ! law is around -0.87. In anoxic zone, it is around -0.35. This
-               ! means a disaggregation constant about 0.5 the value in oxic zones
-               ! -----------------------------------------------------------------
-               zorem  = zremip * trb(ji,jj,jk,jppoc)
-               zofer  = zremip * trb(ji,jj,jk,jpsfe)
-#if ! defined key_kriest
-               zorem2 = zremip * trb(ji,jj,jk,jpgoc)
-               zofer2 = zremip * trb(ji,jj,jk,jpbfe)
-#else
-               zorem2 = zremip * trb(ji,jj,jk,jpnum)
-#endif
-
-               ! Update the appropriate tracers trends
-               ! -------------------------------------
-
-               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer
-#if defined key_kriest
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zorem
-               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) - zorem2
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
-#else
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zorem2 - zorem
-               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) - zorem2
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zofer2 - zofer
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zofer2
-#endif
-
-            END DO
-         END DO
-      END DO
-
-       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
-         WRITE(charout, FMT="('rem3')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
-       ENDIF
-
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zstep   = xstep
-# if defined key_degrad
-               zstep = zstep * facvol(ji,jj,jk)
-# endif
                ! Remineralization rate of BSi depedant on T and saturation
                ! ---------------------------------------------------------
@@ -297 +244 @@
 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
-         WRITE(charout, FMT="('rem4')")
+         WRITE(charout, FMT="('rem3')")
          CALL prt_ctl_trc_info(charout)
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
@@ -315 +262 @@
       END DO
 
-      IF( knt == nrdttrc ) THEN
-          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
-          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
-          !
-          IF( iom_use( "REMIN" ) )  THEN
-              zw3d(:,:,:) = zolimi(:,:,:) * tmask(:,:,:) * zfact !  Remineralisation rate
-              CALL iom_put( "REMIN"  , zw3d )
-          ENDIF
-          IF( iom_use( "DENIT" ) )  THEN
-              zw3d(:,:,:) = denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zfact ! Denitrification
-              CALL iom_put( "DENIT"  , zw3d )
-          ENDIF
-          !
-          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
-       ENDIF
+      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN
+          zrfact2 = 1.e3 * rfact2r
+          CALL iom_put( "REMIN" , zolimi(:,:,:) * tmask(:,:,:) * zrfact2 )  ! Remineralisation rate
+          CALL iom_put( "DENIT" , denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zrfact2  )  ! Denitrification
+      ENDIF
 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
-         WRITE(charout, FMT="('rem6')")
+         WRITE(charout, FMT="('rem4')")
          CALL prt_ctl_trc_info(charout)
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      ztempbac                  )
+      CALL wrk_dealloc( jpi, jpj,      ztempbac )
       CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi )
       !
@@ -357 +294 @@
       !!
       !!----------------------------------------------------------------------
-      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab,   &
+      NAMELIST/nampisrem/ xremik, nitrif, xsirem, xsiremlab, xsilab,   &
       &                   oxymin
       INTEGER :: ios                 ! Local integer output status for namelist read
@@ -374 +311 @@
          WRITE(numout,*) ' Namelist parameters for remineralization, nampisrem'
          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    remineralisation rate of POC              xremip    =', xremip
          WRITE(numout,*) '    remineralization rate of DOC              xremik    =', xremik
          WRITE(numout,*) '    remineralization rate of Si               xsirem    =', xsirem
@@ -386 +322 @@
       denitr  (:,:,:) = 0._wp
       denitnh4(:,:,:) = 0._wp
-      !
+
    END SUBROUTINE p4z_rem_init
 

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

@@ -6 +6 @@
    !! History :   3.5  !  2012-07 (O. Aumont, C. Ethe) Original code
    !!----------------------------------------------------------------------
-#if defined key_pisces
+#if defined key_pisces || defined key_pisces_quota
    !!----------------------------------------------------------------------
    !!   'key_pisces'                                       PISCES bio-model
@@ -42 +42 @@
    REAL(wp), PUBLIC  :: concfediaz  !: Fe half-saturation Cste for diazotrophs 
    REAL(wp)          :: hratio      !: Fe:3He ratio assumed for vent iron supply
+#if defined key_ligand
+   REAL(wp), PUBLIC  :: fep_rats    !: Fep/Fer ratio from sed  sources
+   REAL(wp), PUBLIC  :: fep_rath    !: Fep/Fer ratio from hydro sources
+#endif
 
    LOGICAL , PUBLIC  :: ll_sbc
@@ -72 +76 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdic, rivalk    !: river input fields
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdin, rivdip    !: river input fields
+#if defined key_pisces_quota
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdon, rivdop    !: river input fields
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdoc    !: river input fields
+#endif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: rivdsi    !: river input fields
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:) :: nitdep    !: atmospheric N deposition 
@@ -80 +88 @@
    REAL(wp), PUBLIC :: rivdininput, rivdipinput, rivdsiinput
 
-
    !!* Substitution
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
+   !! $Header:$ 
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -140 +147 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj) 
+                  zcoef = ryyss * cvol(ji,jj,1) 
                   rivalk(ji,jj) =   sf_river(jr_dic)%fnow(ji,jj,1)                                    &
-                     &              * 1.E3        / ( 12. * zcoef + rtrn )
+                     &              * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
+#if defined key_pisces_quota
+                  rivdic(ji,jj) = ( sf_river(jr_dic)%fnow(ji,jj,1) ) &
+                     &              * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
+                  rivdin(ji,jj) = ( sf_river(jr_din)%fnow(ji,jj,1) ) &
+                     &              * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1)
+                  rivdip(ji,jj) = ( sf_river(jr_dip)%fnow(ji,jj,1) ) &
+                     &              * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1)
+                  rivdoc(ji,jj) = ( sf_river(jr_doc)%fnow(ji,jj,1) ) &
+                     &              * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
+                  rivdon(ji,jj) = ( sf_river(jr_don)%fnow(ji,jj,1) ) &
+                     &              * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1)
+                  rivdop(ji,jj) = ( sf_river(jr_dop)%fnow(ji,jj,1) ) &
+                     &              * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1)
+#else
                   rivdic(ji,jj) = ( sf_river(jr_dic)%fnow(ji,jj,1) + sf_river(jr_doc)%fnow(ji,jj,1) ) &
-                     &              * 1.E3         / ( 12. * zcoef + rtrn )
+                     &              * 1.E3 / ( 12. * zcoef + rtrn ) * tmask(ji,jj,1)
                   rivdin(ji,jj) = ( sf_river(jr_din)%fnow(ji,jj,1) + sf_river(jr_don)%fnow(ji,jj,1) ) &
-                     &              * 1.E3 / rno3 / ( 14. * zcoef + rtrn )
+                     &              * 1.E3 / rno3 / ( 14. * zcoef + rtrn ) * tmask(ji,jj,1)
                   rivdip(ji,jj) = ( sf_river(jr_dip)%fnow(ji,jj,1) + sf_river(jr_dop)%fnow(ji,jj,1) ) &
-                     &              * 1.E3 / po4r / ( 31. * zcoef + rtrn )
+                     &              * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1)
+#endif
                   rivdsi(ji,jj) =   sf_river(jr_dsi)%fnow(ji,jj,1)                                    &
-                     &              * 1.E3        / ( 28.1 * zcoef + rtrn )
+                     &              * 1.E3 / ( 28.1 * zcoef + rtrn ) * tmask(ji,jj,1)
                END DO
             END DO
@@ -192 +214 @@
       INTEGER  :: ios                 ! Local integer output status for namelist read
       INTEGER  :: ik50                !  last level where depth less than 50 m
-      INTEGER  :: isrow             ! index for ORCA1 starting row
       REAL(wp) :: zexpide, zdenitide, zmaskt
       REAL(wp) :: ztimes_dust, ztimes_riv, ztimes_ndep 
@@ -208 +229 @@
         &                sn_riverdip, sn_riverdop, sn_riverdsi, sn_ndepo, sn_ironsed, sn_hydrofe, &
         &                ln_dust, ln_solub, ln_river, ln_ndepo, ln_ironsed, ln_ironice, ln_hydrofe,    &
-        &                sedfeinput, dustsolub, icefeinput, wdust, mfrac, nitrfix, diazolight, concfediaz, hratio
+        &                sedfeinput, dustsolub, icefeinput, wdust, mfrac, nitrfix, diazolight, concfediaz,   &
+#if defined key_ligand
+        &                fep_rats, fep_rath,                     &
+#endif
+        &                hratio
       !!----------------------------------------------------------------------
       !
@@ -222 +247 @@
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampissbc in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, nampissbc )
-
-      IF ( ( nn_ice_tr >= 0 ) .AND. ln_ironice ) THEN
-         IF(lwp) THEN
-            WRITE(numout,*) ' ln_ironice incompatible with nn_ice_tr = ', nn_ice_tr
-            WRITE(numout,*) ' Specify your sea ice iron concentration in nampisice instead '
-            WRITE(numout,*) ' ln_ironice is forced to .FALSE. '
-            ln_ironice = .FALSE.
-         ENDIF
-      ENDIF
 
       IF(lwp) THEN
@@ -252 +268 @@
          WRITE(numout,*) '    fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz
          WRITE(numout,*) '    Fe to 3He ratio assumed for vent iron supply hratio  = ', hratio
+#if defined key_ligand
+         WRITE(numout,*) '    Fep/Fer ratio from sed sources                       = ', fep_rats
+         WRITE(numout,*) '    Fep/Fer ratio from sed hydro sources                 = ', fep_rath
+#endif
       END IF
 
@@ -337 +357 @@
          !
          ALLOCATE( rivdic(jpi,jpj), rivalk(jpi,jpj), rivdin(jpi,jpj), rivdip(jpi,jpj), rivdsi(jpi,jpj) ) 
+#if defined key_pisces_quota
+         ALLOCATE( rivdon(jpi,jpj), rivdop(jpi,jpj), rivdoc(jpi,jpj) )
+#endif
          !
          ALLOCATE( sf_river(jpriv), rivinput(jpriv), STAT=ierr1 )           !* allocate and fill sf_river (forcing structure) with sn_river_
@@ -380 +403 @@
          rivalkinput = 0._wp
       END IF 
+
       ! nutrient input from dust
       ! ------------------------
@@ -449 +473 @@
             END DO
          END DO
-         !
+         IF( cp_cfg == 'orca' .AND. jp_cfg == 2 ) THEN
+            ii0 = 176   ;   ii1 =  176        ! Southern Island : Kerguelen
+            ij0 =  37   ;   ij1 =   37  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 119   ;   ii1 =  119        ! South Georgia
+            ij0 =  29   ;   ij1 =   29  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 111   ;   ii1 =  111        ! Falklands
+            ij0 =  35   ;   ij1 =   35  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 168   ;   ii1 =  168        ! Crozet
+            ij0 =  40   ;   ij1 =   40  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 119   ;   ii1 =  119        ! South Orkney
+            ij0 =  28   ;   ij1 =   28  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 140   ;   ii1 =  140        ! Bouvet Island
+            ij0 =  33   ;   ij1 =   33  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 178   ;   ii1 =  178        ! Prince edwards
+            ij0 =  34   ;   ij1 =   34  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 =  43   ;   ii1 =   43        ! Balleny islands
+            ij0 =  21   ;   ij1 =   21  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+         ENDIF
          CALL lbc_lnk( zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
-         !
          DO jk = 1, jpk
             DO jj = 1, jpj
@@ -457 +504 @@
                   zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
                   zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
-                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
+                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 ) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -465 +512 @@
          ironsed(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
-            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday )
+            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday ) 
          END DO
          DEALLOCATE( zcmask)
@@ -483 +530 @@
          CALL iom_close( numhydro )
          !
-         hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp
+         hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp * tmask(:,:,:)
          !
       ENDIF
@@ -519 +566 @@
 
    !!======================================================================
-END MODULE p4zsbc
+END MODULE  p4zsbc

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

@@ -75 +75 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsbio3, zwsbio4, zwscal
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zirondep, zsoufer
+#if defined key_ligand
+      REAL(wp) ::  zwssfep
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsfep
+#endif
       !!---------------------------------------------------------------------
       !
@@ -85 +89 @@
       CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal )
       CALL wrk_alloc( jpi, jpj, jpk, zsoufer )
+#if defined key_ligand
+      CALL wrk_alloc( jpi, jpj, zwsfep )
+#endif
 
       zdenit2d(:,:) = 0.e0
@@ -186 +193 @@
       IF( ln_ironsed ) THEN
          tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
+#if defined key_ligand
+         trn(:,:,:,jpfep) = trn(:,:,:,jpfep) + (ironsed(:,:,:) * fep_rats ) * rfact2
+#endif
          !
          IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
@@ -195 +205 @@
       IF( ln_hydrofe ) THEN
          tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
+#if defined key_ligand
+         trn(:,:,:,jpfep) = trn(:,:,:,jpfep) + ( hydrofe(:,:,:) * fep_rath ) * rfact2
+         trn(:,:,:,jplgw) = trn(:,:,:,jplgw) + ( hydrofe(:,:,:) * 0.5 ) * rfact2
+#endif
          !
          IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) )   &
@@ -210 +224 @@
             zwscal (ji,jj) = MIN( 0.99 * zdep, wscal (ji,jj,ikt) )
             zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
+#if defined key_ligand
+            zwsfep(ji,jj)  = MIN( 0.99 * zdep, wsfep(ji,jj,ikt)  )
+#endif
          END DO
       END DO
@@ -308 +325 @@
             zws4 = zwsbio4(ji,jj) * zdep
             zws3 = zwsbio3(ji,jj) * zdep
+#if defined key_ligand
+            zwssfep = zwsfep(ji,jj) * zdep
+#endif
             zrivno3 = 1. - zbureff(ji,jj)
 # if ! defined key_kriest
@@ -315 +335 @@
             tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3
             zwstpoc              = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3
+#   if defined key_ligand
+            trn(ji,jj,ikt,jpfep) = trn(ji,jj,ikt,jpfep) - trn(ji,jj,ikt,jpfep) * zwssfep
+#   endif
 # else
             tra(ji,jj,ikt,jpnum) = tra(ji,jj,ikt,jpnum) - trb(ji,jj,ikt,jpnum) * zws4 
@@ -320 +343 @@
             tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3
             zwstpoc = trb(ji,jj,ikt,jppoc) * zws3 
+#   if defined key_ligand
+            trn(ji,jj,ikt,jpfep) = trn(ji,jj,ikt,jpfep) - trn(ji,jj,ikt,jpfep) * zwssfep
+#   endif
 # endif
 
@@ -407 +433 @@
       CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal )
       CALL wrk_dealloc( jpi, jpj, jpk, zsoufer )
+#if defined key_ligand
+      CALL wrk_dealloc( jpi, jpj, zwsfep )
+#endif
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_sed')

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

@@ -40 +40 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer2           !: Big iron sinking fluxes
 #endif
+#if defined key_ligand
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfep      !: Fep sinking fluxes
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsfep
+#endif
 
    INTEGER  :: ik100
@@ -91 +95 @@
       INTEGER  ::   ji, jj, jk, jit
       INTEGER  ::   iiter1, iiter2
-      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
-      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
       REAL(wp) ::   zfact, zwsmax, zmax, zstep
       CHARACTER (len=25) :: charout
@@ -100 +102 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_sink')
+
       !
       !    Sinking speeds of detritus is increased with depth as shown
@@ -108 +111 @@
             DO ji = 1,jpi
                zmax  = MAX( heup(ji,jj), hmld(ji,jj) )
-               zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / 5000._wp
-               wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
+               zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / wsbio2scale
+               wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
             END DO
          END DO
@@ -117 +120 @@
       wsbio3(:,:,:) = wsbio
       wscal (:,:,:) = wsbio4(:,:,:)
+#if defined key_ligand
+      wsfep (:,:,:) = wfep
+#endif
       !
       ! OA This is (I hope) a temporary solution for the problem that may 
@@ -159 +165 @@
                  wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
                  wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) )
+#if defined key_ligand
+                 wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
+#endif
                ENDIF
             END DO
@@ -172 +181 @@
       sinksil (:,:,:) = 0.e0
       sinkfer2(:,:,:) = 0.e0
+#if defined key_ligand
+      sinkfep(:,:,:) = 0.e0
+#endif
 
       !   Compute the sedimentation term using p4zsink2 for all the sinking particles
@@ -178 +190 @@
         CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 )
         CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 )
+#if defined key_ligand
+        CALL p4z_sink2( wsfep , sinkfep , jpfep, iiter1 )
+#endif
       END DO
 
@@ -186 +201 @@
         CALL p4z_sink2( wscal , sinkcal , jpcal, iiter2 )
       END DO
-
-      !  Exchange between organic matter compartments due to coagulation/disaggregation
-      !  ---------------------------------------------------
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               zstep = xstep 
-# if defined key_degrad
-               zstep = zstep * facvol(ji,jj,jk)
-# endif
-               zfact = zstep * xdiss(ji,jj,jk)
-               !  Part I : Coagulation dependent on turbulence
-               zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
-               zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
-
-               ! Part II : Differential settling
-
-               !  Aggregation of small into large particles
-               zagg3 =  47.1 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
-               zagg4 =  3.3  * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
-
-               zagg   = zagg1 + zagg2 + zagg3 + zagg4
-               zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn )
-
-               ! Aggregation of DOC to POC : 
-               ! 1st term is shear aggregation of DOC-DOC
-               ! 2nd term is shear aggregation of DOC-POC
-               ! 3rd term is differential settling of DOC-POC
-               zaggdoc  = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
-               &            + 2.4 * zstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
-               ! transfer of DOC to GOC : 
-               ! 1st term is shear aggregation
-               ! 2nd term is differential settling 
-               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
-               ! tranfer of DOC to POC due to brownian motion
-               zaggdoc3 =  ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)
-
-               !  Update the trends
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc + zaggdoc3
-               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
-               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
-               !
-            END DO
-         END DO
-      END DO
-
 
      ! Total carbon export per year
@@ -341 +307 @@
       !
       INTEGER  :: ji, jj, jk, jit, niter1, niter2
-      REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggsi, zaggsh
-      REAL(wp) :: zagg , zaggdoc, zaggdoc1, znumdoc
       REAL(wp) :: znum , zeps, zfm, zgm, zsm
       REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
-      REAL(wp) :: zval1, zval2, zval3, zval4
+      REAL(wp) :: zval1, zval2, zval3
       REAL(wp) :: zfact
       INTEGER  :: ik1
@@ -395 +359 @@
 
       wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp )
+#if defined key_ligand
+      wsfep (:,:,:) = wfep
+#endif
 
       !   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
@@ -404 +371 @@
       sinkfer (:,:,:) = 0.e0
       sinksil (:,:,:) = 0.e0
+#if defined key_ligand
+      sinkfep(:,:,:) = 0.e0
+#endif
 
      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
@@ -416 +386 @@
         CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )
         CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )
+#if defined key_ligand
+        CALL p4z_sink2( wsfep,  sinkfep , jpfep, iiter1 )
+#endif
       END DO
 
@@ -422 +395 @@
       END DO
 
-     !  Exchange between organic matter compartments due to coagulation/disaggregation
-     !  ---------------------------------------------------
-
-      zval1 = 1. + xkr_zeta
-      zval2 = 1. + xkr_eta
-      zval3 = 3. + xkr_eta
-      zval4 = 4. + xkr_eta
-
-      DO jk = 1,jpkm1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
-
-                  znum = trb(ji,jj,jk,jppoc)/(trb(ji,jj,jk,jpnum)+rtrn) / xkr_massp
-                  !-------------- To avoid sinking speed over 50 m/day -------
-                  znum  = min(xnumm(jk),znum)
-                  znum  = MAX( 1.1,znum)
-                  !------------------------------------------------------------
-                  zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
-                  zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
-                  zdiv1 = MAX( 1.e-4, ABS( zeps - 4.   ) ) * SIGN( 1., zeps - 4.    )
-                  zdiv2 = zeps - 2.
-                  zdiv3 = zeps - 3.
-                  zdiv4 = zeps - zval2
-                  zdiv5 = 2.* zeps - zval4
-                  zfm   = xkr_frac**( 1.- zeps )
-                  zsm   = xkr_frac**xkr_eta
-
-                  !    Part I : Coagulation dependant on turbulence
-                  !    ----------------------------------------------
-
-                  zagg1 =  0.163 * trb(ji,jj,jk,jpnum)**2               &
-                     &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    &
-                     &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    &
-                     &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
-                     &            * (zeps-1.)**2/(zdiv2*zdiv3)) 
-                  zagg2 =  2*0.163*trb(ji,jj,jk,jpnum)**2*zfm*                       &
-                     &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2          &
-                     &                    *xkr_mass_min*(zeps-1.)/zdiv2                 &
-                     &                    +xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3)    &
-                     &                    +xkr_mass_min**3*(zeps-1)/zdiv1)                  &
-                     &                    -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/           &
-                     &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))    
-
-                  zagg3 =  0.163*trb(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3  
-                  
-                 !    Aggregation of small into large particles
-                 !    Part II : Differential settling
-                 !    ----------------------------------------------
-
-                  zagg4 =  2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2*                       &
-                     &                 xkr_wsbio_min*(zeps-1.)**2                         &
-                     &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      &
-                     &                 -(1.-zfm)/(zdiv*(zeps-1.)))-                       &
-                     &                 ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2)     &
-                     &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )   
-
-                  zagg5 =   2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2                         &
-                     &                 *(zeps-1.)*zfm*xkr_wsbio_min                        &
-                     &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         &
-                     &                 /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2)    &
-                     &                 /zdiv)  
-
-                  !
-                  !     Fractionnation by swimming organisms
-                  !     ------------------------------------
-
-                  zfract = 2.*3.141*0.125*trb(ji,jj,jk,jpmes)*12./0.12/0.06**3*trb(ji,jj,jk,jpnum)  &
-                    &      * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2  &
-                    &      * 10000.*xstep
-
-                  !     Aggregation of DOC to small particles
-                  !     --------------------------------------
-
-                  zaggdoc = 0.83 * trb(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)   &
-                     &        + 0.005 * 231. * trb(ji,jj,jk,jpdoc) * xstep * trb(ji,jj,jk,jpdoc)
-                  zaggdoc1 = 271. * trb(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)  &
-                     &  + 0.02 * 16706. * trb(ji,jj,jk,jppoc) * xstep * trb(ji,jj,jk,jpdoc)
-
-# if defined key_degrad
-                   zagg1   = zagg1   * facvol(ji,jj,jk)                 
-                   zagg2   = zagg2   * facvol(ji,jj,jk)                 
-                   zagg3   = zagg3   * facvol(ji,jj,jk)                 
-                   zagg4   = zagg4   * facvol(ji,jj,jk)                 
-                   zagg5   = zagg5   * facvol(ji,jj,jk)                 
-                   zaggdoc = zaggdoc * facvol(ji,jj,jk)                 
-                   zaggdoc1 = zaggdoc1 * facvol(ji,jj,jk)
-# endif
-                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
-                  zaggsi = ( zagg4 + zagg5 ) * xstep / 10.
-                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
-                  !
-                  znumdoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )
-                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1
-                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg
-                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc1
-
-               ENDIF
-            END DO
-         END DO
-      END DO
-
-     ! Total primary production per year
-     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) )
+     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
+        &   t_oce_co2_exp = glob_sum( sinking(:,:,ik100)  * e1e2t(:,:) * tmask(:,:,1) )
      !
      IF( lk_iomput ) THEN
@@ -755 +626 @@
       ik100 = 10        !  last level where depth less than 100 m
       DO jk = jpkm1, 1, -1
-         IF( gdept_1d(jk) > 100. )  iksed = jk - 1
+         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1
       END DO
       IF (lwp) WRITE(numout,*)
@@ -897 +768 @@
          &      sinkfer2(jpi,jpj,jpk)                                             ,     &                
 #endif
+#if defined key_ligand
+         &      wsfep(jpi,jpj,jpk)   , sinkfep(jpi,jpj,jpk)                         ,     &
+#endif
          &      sinkfer(jpi,jpj,jpk)                                              , STAT=p4z_sink_alloc )                
          !

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

@@ -19 +19 @@
    USE p4zbio          !  Biological model
    USE p4zche          !  Chemical model
-   USE p4zlys          !  Calcite saturation
    USE p4zflx          !  Gas exchange
    USE p4zsbc          !  External source of nutrients
@@ -85 +84 @@
         CALL p4z_che                              ! initialize the chemical constants
         !
-        IF( .NOT. ln_rsttr ) THEN  ;   CALL p4z_ph_ini   !  set PH at kt=nit000 
+        IF( .NOT. ln_rsttr ) THEN  ;   CALL ahini_for_at(hi)   !  set PH at kt=nit000
         ELSE                       ;   CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields 
         ENDIF
@@ -133 +132 @@
          !
          CALL p4z_bio( kt, jnt )   ! Biology
+         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
          CALL p4z_sed( kt, jnt )   ! Sedimentation
-         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
-         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
          !
          xnegtr(:,:,:) = 1.e0
@@ -216 +214 @@
       !!                       natkriest ("key_kriest")
       !!----------------------------------------------------------------------
-      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, niter1max, niter2max
+#if defined key_ligand
+      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    &
+      &                   niter1max, niter2max, wfep, ldocp, ldocz, lthet
+#else
+      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale, niter1max, niter2max
+#endif
 #if defined key_kriest
       NAMELIST/nampiskrp/ xkr_eta, xkr_zeta, xkr_ncontent, xkr_mass_min, xkr_mass_max
@@ -241 +244 @@
          WRITE(numout,*) '    Fe/C in zooplankton                       ferat3    =', ferat3
          WRITE(numout,*) '    Big particles sinking speed               wsbio2    =', wsbio2
+         WRITE(numout,*) '    Big particles maximum sinking speed       wsbio2max =', wsbio2max
+         WRITE(numout,*) '    Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
          WRITE(numout,*) '    Maximum number of iterations for POC      niter1max =', niter1max
          WRITE(numout,*) '    Maximum number of iterations for GOC      niter2max =', niter2max
+#if defined key_ligand
+         WRITE(numout,*) '    FeP sinking speed                             wfep  =', wfep
+         WRITE(numout,*) '    Phyto ligand production per unit doc          ldocp =', ldocp
+         WRITE(numout,*) '    Zoo ligand production per unit doc            ldocz =', ldocz
+         WRITE(numout,*) '    Proportional loss of ligands due to Fe uptake lthet =', lthet
+#endif
       ENDIF
 
@@ -310 +321 @@
    END SUBROUTINE p4z_sms_init
 
-   SUBROUTINE p4z_ph_ini
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE p4z_ini_ph  ***
-      !!
-      !!  ** Purpose : Initialization of chemical variables of the carbon cycle
-      !!---------------------------------------------------------------------
-      INTEGER  ::  ji, jj, jk
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
-      !!---------------------------------------------------------------------
-
-      ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???)
-      ! --------------------------------------------------------
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ztmas   = tmask(ji,jj,jk)
-               ztmas1  = 1. - tmask(ji,jj,jk)
-               zcaralk = trb(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
-               zco3    = ( zcaralk - trb(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
-               zbicarb = ( 2. * trb(ji,jj,jk,jpdic) - zcaralk )
-               hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
-            END DO
-         END DO
-     END DO
-     !
-   END SUBROUTINE p4z_ph_ini
-
    SUBROUTINE p4z_rst( kt, cdrw )
       !!---------------------------------------------------------------------
@@ -350 +333 @@
       INTEGER         , INTENT(in) ::   kt         ! ocean time-step
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
-      !
-      INTEGER  ::  ji, jj, jk
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
       !!---------------------------------------------------------------------
 
@@ -365 +344 @@
             CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
          ELSE
-!            hi(:,:,:) = 1.e-9 
-            CALL p4z_ph_ini
+            CALL ahini_for_at(hi)
          ENDIF
          CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
@@ -550 +528 @@
             &                    + trn(:,:,:,jpbfe)                     &
 #endif
+#if defined key_ligand
+            &                    + trn(:,:,:,jpfep)                     &
+#endif
             &                    + trn(:,:,:,jpsfe)                     &
             &                    + trn(:,:,:,jpzoo) * ferat3            &

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/par_pisces.F90

@@ -18 +18 @@
    !!---------------------------------------------------------------------
    LOGICAL, PUBLIC, PARAMETER ::   lk_pisces     = .TRUE.  !: PISCES flag 
-   LOGICAL, PUBLIC, PARAMETER ::   lk_p4z        = .FALSE. !: p4z flag 
+   INTEGER, PUBLIC, PARAMETER ::   nn_p4z        =  1      !: p4z flag 
    INTEGER, PUBLIC, PARAMETER ::   jp_pisces     =  6      !: number of passive tracers
    INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  =  19     !: additional 2d output 
@@ -36 +36 @@
    INTEGER, PUBLIC, PARAMETER ::   jpkbm1    = jpkb - 1  !: first vertical layers where biology is active
 
-#elif defined key_pisces  &&  defined key_kriest
+#elif defined key_pisces
    !!---------------------------------------------------------------------
    !!   'key_pisces' & 'key_kriest'                 PISCES bio-model + ???
    !!---------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_pisces     = .TRUE.  !: PISCES flag 
-   LOGICAL, PUBLIC, PARAMETER ::   lk_p4z        = .TRUE. !: p4z flag 
-   LOGICAL, PUBLIC, PARAMETER ::   lk_kriest     = .TRUE.  !: Kriest flag 
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces     =  23     !: number of passive tracers
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  =  13     !: additional 2d output 
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  =  18     !: additional 3d output 
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd =   1     !: number of sms trends for PISCES
+   LOGICAL, PUBLIC, PARAMETER ::   lk_pisces      = .TRUE.  !: PISCES flag 
+   INTEGER, PUBLIC, PARAMETER ::   nn_p4z         =  2      !: p4z flag 
+#  if defined key_kriest
+   LOGICAL, PUBLIC, PARAMETER ::   lk_kriest      = .TRUE.  !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_kriest      =  1 
+   INTEGER, PUBLIC, PARAMETER ::   jp_kriest_diag =  7
+#  else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_kriest      = .FALSE.  !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_kriest      =  2 
+   INTEGER, PUBLIC, PARAMETER ::   jp_kriest_diag =  0
+#  endif
+#  if defined key_ligand
+   LOGICAL, PUBLIC, PARAMETER ::   lk_ligand      = .TRUE.  !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_ligand      =  2 
+   INTEGER, PUBLIC, PARAMETER ::   jp_ligand_diag =  0
+#  else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_ligand      = .FALSE.  !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_ligand      =  0     
+   INTEGER, PUBLIC, PARAMETER ::   jp_ligand_diag =  0
+#  endif
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces      =  22 + jp_kriest + jp_ligand            !: number of passive tracers
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d   =  13                                    !: additional 2d output 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d   =  13 + jp_kriest_diag + jp_ligand_diag  !: additional 3d output 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd  =   1                                    !: number of sms trends for PISCES
 
    ! assign an index in trc arrays for each LOBSTER prognostic variables
@@ -63 +80 @@
    INTEGER, PUBLIC, PARAMETER ::   jpdia = 11    !: Diatoms Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpmes = 12    !: Mesozooplankton Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpdsi = 13    !: Diatoms Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdsi = 13    !: (big) Silicate Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpfer = 14    !: Iron Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpnum = 15    !: Big iron particles Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpsfe = 16    !: number of particulate organic phosphate concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpdfe = 17    !: Diatoms iron Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpgsi = 18    !: (big) Silicate Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpnfe = 19    !: Nano iron Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpnch = 20    !: Nano Chlorophyll Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpdch = 21    !: Diatoms Chlorophyll Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpno3 = 22    !: Nitrates Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpnh4 = 23    !: Ammonium Concentration
-
-#elif defined key_pisces
-   !!---------------------------------------------------------------------
-   !!   'key_pisces'   :                         standard PISCES bio-model
+   INTEGER, PUBLIC, PARAMETER ::   jpsfe = 15    !: number of particulate organic phosphate concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdfe = 16    !: Diatoms iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpgsi = 17    !: Diatoms Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnfe = 18    !: Nano iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnch = 19    !: Nano Chlorophyll Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdch = 20    !: Diatoms Chlorophyll Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpno3 = 21    !: Nitrates Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnh4 = 22    !: Ammonium Concentration
+#  if defined key_kriest
+   INTEGER, PUBLIC, PARAMETER ::   jpnum = 23    !: Number of particles in the aggregates
+#  else
+   INTEGER, PUBLIC, PARAMETER ::   jpgoc = 23    !: Big OM particles Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpbfe = 24    !: Big iron particles Concentration
+#  endif
+#  if defined key_ligand
+   INTEGER, PUBLIC, PARAMETER ::   jplgw = 22 + jp_kriest + 1    !: Weak Ligands
+   INTEGER, PUBLIC, PARAMETER ::   jpfep = 22 + jp_kriest + 2    !: Fe nanoparticle
+#  endif
+
+#elif defined key_pisces_quota 
+   !!---------------------------------------------------------------------
+   !!   'key_pisces' & 'key_kriest'                 PISCES bio-model + ???
    !!---------------------------------------------------------------------
    LOGICAL, PUBLIC, PARAMETER ::   lk_pisces     = .TRUE.  !: PISCES flag 
-   LOGICAL, PUBLIC, PARAMETER ::   lk_p4z        = .TRUE.  !: p4z flag 
+   INTEGER, PUBLIC, PARAMETER ::   nn_p4z        =  3  !: p4z flag 
+#  if defined key_kriest
+   LOGICAL, PUBLIC, PARAMETER ::   lk_kriest     = .TRUE. !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_kriest      =  1
+   INTEGER, PUBLIC, PARAMETER ::   jp_kriest_diag =  7
+#  else
    LOGICAL, PUBLIC, PARAMETER ::   lk_kriest     = .FALSE. !: Kriest flag 
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces     = 24      !: number of PISCES passive tracers
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  = 13      !: additional 2d output 
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d  = 11      !: additional 3d output 
-   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd =  1      !: number of sms trends for PISCES
+   INTEGER, PUBLIC, PARAMETER ::   jp_kriest      =  4
+   INTEGER, PUBLIC, PARAMETER ::   jp_kriest_diag =  7
+#  endif
+#  if defined key_ligand
+   LOGICAL, PUBLIC, PARAMETER ::   lk_ligand      = .TRUE.  !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_ligand      =  2
+   INTEGER, PUBLIC, PARAMETER ::   jp_ligand_diag =  0
+#  else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_ligand      = .FALSE.  !: Kriest flag 
+   INTEGER, PUBLIC, PARAMETER ::   jp_ligand      =  0
+   INTEGER, PUBLIC, PARAMETER ::   jp_ligand_diag =  0
+#  endif
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces      =  35 + jp_kriest + jp_ligand            !: number of passive tracers
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d   =  13     !: additional 2d output 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_3d   =  13 + jp_kriest_diag + jp_ligand_diag  !: additional 3d output 
+   INTEGER, PUBLIC, PARAMETER ::   jp_pisces_trd  =   1     !: number of sms trends for PISCES
 
    ! assign an index in trc arrays for each LOBSTER prognostic variables
@@ -102 +145 @@
    INTEGER, PUBLIC, PARAMETER ::   jpdia = 11    !: Diatoms Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpmes = 12    !: Mesozooplankton Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpdsi = 13    !: Diatoms Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdsi = 13    !: (big) Silicate Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpfer = 14    !: Iron Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpbfe = 15    !: Big iron particles Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpgoc = 16    !: big particulate organic phosphate concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpsfe = 17    !: Small iron particles Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpdfe = 18    !: Diatoms iron Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpgsi = 19    !: (big) Silicate Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpnfe = 20    !: Nano iron Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpnch = 21    !: Nano Chlorophyll Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpdch = 22    !: Diatoms Chlorophyll Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpno3 = 23    !: Nitrates Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpnh4 = 24    !: Ammonium Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpsfe = 15    !: number of particulate organic phosphate concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdfe = 16    !: Diatoms iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpgsi = 17    !: Diatoms Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnfe = 18    !: Nano iron Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnch = 19    !: Nano Chlorophyll Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdch = 20    !: Diatoms Chlorophyll Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpno3 = 21    !: Nitrates Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnh4 = 22    !: Ammonium Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdon = 23    !: dissolved organic nitrogen concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdop = 24    !: dissolved organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jppon = 25    !: small particulate organic nitrogen concentration
+   INTEGER, PUBLIC, PARAMETER ::   jppop = 26    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnph = 27    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jppph = 28    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpndi = 29    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jppdi = 30    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jppic = 31    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpnpi = 32    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpppi = 33    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jppfe = 34    !: small particulate organic phosphorus concentration
+   INTEGER, PUBLIC, PARAMETER ::   jppch = 35    !: small particulate organic phosphorus concentration
+#  if defined key_kriest
+   INTEGER, PUBLIC, PARAMETER ::   jpnum = 36    !: Number of particles in the aggregates
+#  else
+   INTEGER, PUBLIC, PARAMETER ::   jpgoc = 36    !: Big carbon particles Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpgon = 37    !: Big nitrogen particles Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpgop = 38    !: Big phosphorus particles Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpbfe = 39    !: Big iron particles Concentration
+#  endif
+#  if defined key_ligand
+   INTEGER, PUBLIC, PARAMETER ::   jplgw = 35 + jp_kriest + 1    !: Weak Ligands
+   INTEGER, PUBLIC, PARAMETER ::   jpfep = 35 + jp_kriest + 2    !: Fe nanoparticle
+#  endif
 
 #else
@@ -120 +186 @@
    !!---------------------------------------------------------------------
    LOGICAL, PUBLIC, PARAMETER ::   lk_pisces     = .FALSE.  !: PISCES flag 
-   LOGICAL, PUBLIC, PARAMETER ::   lk_p4z        = .FALSE.  !: p4z flag 
+   INTEGER, PUBLIC, PARAMETER ::   nn_p4z        =  0  !: p4z flag 
    INTEGER, PUBLIC, PARAMETER ::   jp_pisces     =  0       !: No CFC tracers
    INTEGER, PUBLIC, PARAMETER ::   jp_pisces_2d  =  0       !: No CFC additional 2d output arrays 

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90

@@ -6 +6 @@
    !! History :   1.0  !  2000-02 (O. Aumont) original code
    !!             3.2  !  2009-04 (C. Ethe & NEMO team) style
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
-#if defined key_pisces || defined key_pisces_reduced 
+#if defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota 
    !!----------------------------------------------------------------------
-   !!   'key_pisces'                                         PISCES model
+   !!   'key_pisces*'                                         PISCES model
    !!----------------------------------------------------------------------
    USE par_oce
@@ -22 +23 @@
 
    !!*  Biological fluxes for light : variables shared by pisces & lobster
-   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  neln  !: number of T-levels + 1 in the euphotic layer
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  heup  !: euphotic layer depth
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  etot  !: par (photosynthetic available radiation)
-   !
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  xksi  !:  LOBSTER : zooplakton closure
-   !                                                       !:  PISCES  : silicon dependant half saturation
+   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  neln    !: number of T-levels + 1 in the euphotic layer
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  heup    !: euphotic layer depth
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  heup_01 !: Absolute euphotic layer depth
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  etot    !: par (photosynthetic available radiation)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  xksi    !:  LOBSTER : zooplakton closure
+   !                                                         !:  PISCES  : silicon dependant half saturation
 
-#if defined key_pisces 
+#if defined key_pisces || defined key_pisces_quota 
    !!*  Time variables
    INTEGER  ::   nrdttrc           !: ???
@@ -38 +39 @@
    REAL(wp) ::   ryyss             !: number of seconds per year 
    REAL(wp) ::   r1_ryyss          !: inverse number of seconds per year 
-
 
    !!*  Biological parameters 
@@ -49 +49 @@
    REAL(wp) ::   o2nit             !: ???
    REAL(wp) ::   wsbio, wsbio2     !: ???
+   REAL(wp) ::   wsbio2max         !: ???
+   REAL(wp) ::   wsbio2scale       !: ???
    REAL(wp) ::   xkmort            !: ???
    REAL(wp) ::   ferat3            !: ???
+#if defined key_ligand
+   REAL(wp) ::   wfep              !: ???
+   REAL(wp) ::   ldocp           !: ???
+   REAL(wp) ::   ldocz           !: ???
+   REAL(wp) ::   lthet           !: ???
+#endif
+
 
    !!*  diagnostic parameters 
@@ -68 +77 @@
    !!*  Biological fluxes for primary production
    REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:)  ::   xksimax    !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanono3   !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatno3   !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanonh4   !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatnh4   !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanopo4   !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatpo4   !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimphy    !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdia    !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concdfe    !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concnfe    !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimnfe    !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdfe    !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimsi     !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   biron      !: bioavailable fraction of iron
-
+#if defined key_ligand
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   plig       !: proportion of iron organically complexed
+#endif
 
    !!*  SMS for the organic matter
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xfracal    !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nitrfac    !: ??
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xlimbac    !: ??
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xlimbacl   !: ??
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   orem       !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xdiss      !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodcal    !: Calcite production
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodpoc    !: Calcite production
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   conspoc    !: Calcite production
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodgoc    !: Calcite production
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   consgoc    !: Calcite production
+
 
    !!* Variable for chemistry of the CO2 cycle
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak13       !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak23       !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksp       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akw3       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   borat      !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi         !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   excess     !: ???
@@ -116 +115 @@
 
 #endif
+
+# if defined key_pisces_quota
+   !!*  Biological parameters 
+   REAL(wp) ::   no3rat3           !: ???
+   REAL(wp) ::   po4rat3           !: ???
+
+   !!*  SMS for the organic matter
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizen      !: size of diatoms 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizep      !: size of diatoms 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sized      !: size of diatoms 
+
+#endif
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -128 +139 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_warn
+#if defined key_pisces_quota
+      INTEGER ::   ierr(6)        ! Local variables
+#else
       INTEGER ::   ierr(5)        ! Local variables
+#endif
       !!----------------------------------------------------------------------
       ierr(:) = 0
       !*  Biological fluxes for light : shared variables for pisces & lobster
-      ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj), xksi(jpi,jpj), STAT=ierr(1) )
+      ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj),    &
+      &         heup_01(jpi,jpj), xksi(jpi,jpj), STAT=ierr(1) )
       !
-#if defined key_pisces
+#if defined key_pisces || defined key_pisces_quota
       !*  Biological fluxes for primary production
       ALLOCATE( xksimax(jpi,jpj)     , biron   (jpi,jpj,jpk),       &
-         &      xnanono3(jpi,jpj,jpk), xdiatno3(jpi,jpj,jpk),       &
-         &      xnanonh4(jpi,jpj,jpk), xdiatnh4(jpi,jpj,jpk),       &
-         &      xnanopo4(jpi,jpj,jpk), xdiatpo4(jpi,jpj,jpk),       &
-         &      xlimphy (jpi,jpj,jpk), xlimdia (jpi,jpj,jpk),       &
-         &      xlimnfe (jpi,jpj,jpk), xlimdfe (jpi,jpj,jpk),       &
-         &      xlimsi  (jpi,jpj,jpk), concdfe (jpi,jpj,jpk),       &
-         &      concnfe (jpi,jpj,jpk),                           STAT=ierr(2) ) 
+#if defined key_ligand
+      &         plig(jpi,jpj,jpk)    ,                              &
+#endif
+         &      STAT=ierr(2) )
          !
       !*  SMS for the organic matter
-      ALLOCATE( xfracal (jpi,jpj,jpk), nitrfac(jpi,jpj,jpk),       &
-         &      xlimbac (jpi,jpj,jpk), xdiss  (jpi,jpj,jpk),       & 
-         &      xlimbacl(jpi,jpj,jpk), prodcal(jpi,jpj,jpk),     STAT=ierr(3) )
+      ALLOCATE( xfracal (jpi,jpj,jpk), nitrfac (jpi,jpj,jpk),      &
+         &      orem    (jpi,jpj,jpk),                             &
+         &      prodcal(jpi,jpj,jpk),  xdiss   (jpi,jpj,jpk),      & 
+         &      prodpoc(jpi,jpj,jpk) , conspoc(jpi,jpj,jpk),       &
+         &      prodgoc(jpi,jpj,jpk) , consgoc(jpi,jpj,jpk),     STAT=ierr(3) )
 
       !* Variable for chemistry of the CO2 cycle
-      ALLOCATE( akb3(jpi,jpj,jpk)    , ak13  (jpi,jpj,jpk) ,       &
+      ALLOCATE( ak13  (jpi,jpj,jpk) ,                              &
          &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,       &
-         &      akw3(jpi,jpj,jpk)    , borat (jpi,jpj,jpk) ,       &
          &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,     STAT=ierr(4) )
          !
@@ -160 +174 @@
          !
 #endif
+#if defined key_pisces_quota
+      !*  Size of phytoplankton cells
+      ALLOCATE( sizen   (jpi,jpj,jpk), sizep   (jpi,jpj,jpk),      &
+         &      sized   (jpi,jpj,jpk), STAT=ierr(6) )
+         !
+#endif
+
       !
       sms_pisces_alloc = MAXVAL( ierr )

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcice_pisces.F90

@@ -37 +37 @@
       !!----------------------------------------------------------------------
 
-      IF( lk_p4z ) THEN  ;   CALL p4z_ice_ini   !  PISCES
-      ELSE               ;   CALL p2z_ice_ini   !  LOBSTER
+      IF( nn_p4z == 1 ) THEN  ;   CALL p2z_ice_ini   !  PISCES
+      ELSE                    ;   CALL p4z_ice_ini   !  LOBSTER
       ENDIF
 
@@ -128 +128 @@
       zpisc(jpno3,1) =  5.79e-6_wp / rno3 
       zpisc(jpnh4,1) =  3.22e-7_wp / rno3
+#if defined key_pisces_quota
+      zpisc(jppic,1) =  9.57e-8_wp
+      zpisc(jpnpi,1) =  9.57e-8_wp
+      zpisc(jpppi,1) =  9.57e-8_wp
+      zpisc(jppfe,1) =  1.76e-11_wp
+      zpisc(jppch,1) =  1.67e-7_wp
+      zpisc(jpnph,1) =  9.57e-8_wp
+      zpisc(jppph,1) =  9.57e-8_wp
+      zpisc(jpndi,1) =  4.24e-7_wp
+      zpisc(jppdi,1) =  4.24e-7_wp
+      zpisc(jppon,1) =  9.57e-8_wp
+      zpisc(jppop,1) =  9.57e-8_wp
+      zpisc(jpdon,1) =  2.04e-5_wp
+      zpisc(jpdop,1) =  2.04e-5_wp
+      zpisc(jpgon,1) =  5.23e-8_wp
+      zpisc(jpgop,1) =  5.23e-8_wp
+#endif
 
       !--- Arctic specificities (dissolved inorganic & DOM)
@@ -158 +175 @@
       zpisc(jpno3,2) =  3.51e-06_wp / rno3 
       zpisc(jpnh4,2) =  6.15e-08_wp / rno3 
+#if defined key_pisces_quota
+      zpisc(jppic,1) =  5.25e-7_wp
+      zpisc(jpnpi,1) =  5.25e-7_wp
+      zpisc(jpppi,1) =  5.25e-7_wp
+      zpisc(jppfe,1) =  1.75e-11_wp
+      zpisc(jppch,1) =  1.46e-07_wp
+      zpisc(jpnph,1) =  5.25e-7_wp
+      zpisc(jppph,1) =  5.25e-7_wp
+      zpisc(jpndi,1) =  7.75e-7_wp
+      zpisc(jppdi,1) =  7.75e-7_wp
+      zpisc(jppon,1) =  4.05e-7_wp
+      zpisc(jppop,1) =  4.05e-7_wp
+      zpisc(jpdon,1) =  6.00e-6_wp
+      zpisc(jpdop,1) =  6.00e-6_wp
+      zpisc(jpgon,1) =  2.84e-8_wp
+      zpisc(jpgop,1) =  2.84e-8_wp
+#endif
+
 
       !--- Antarctic specificities (dissolved inorganic & DOM)
@@ -188 +223 @@
       zpisc(jpno3,3) =  2.64e-5_wp / rno3  
       zpisc(jpnh4,3) =  3.39e-7_wp / rno3  
+#if defined key_pisces_quota
+      zpisc(jppic,1) =  8.10e-7_wp
+      zpisc(jpnpi,1) =  8.10e-7_wp
+      zpisc(jpppi,1) =  8.10e-7_wp 
+      zpisc(jppfe,1) =  1.48e-11_wp
+      zpisc(jppch,1) =  2.02e-7_wp
+      zpisc(jpnph,1) =  9.57e-8_wp
+      zpisc(jppph,1) =  9.57e-8_wp
+      zpisc(jpndi,1) =  5.77e-7_wp
+      zpisc(jppdi,1) =  5.77e-7_wp
+      zpisc(jppon,1) =  1.13e-6_wp
+      zpisc(jppop,1) =  1.13e-6_wp
+      zpisc(jpdon,1) =  7.02e-6_wp
+      zpisc(jpdop,1) =  7.02e-6_wp
+      zpisc(jpgon,1) =  2.89e-8_wp
+      zpisc(jpgop,1) =  2.89e-8_wp
+#endif
+
 
       !--- Baltic Sea particular case for ORCA configurations
@@ -218 +271 @@
       zpisc(jpno3,4) = 5.36e-5_wp / rno3
       zpisc(jpnh4,4) = 7.18e-7_wp / rno3
+#if defined key_pisces_quota
+      zpisc(jppic,1) =  6.64e-7_wp
+      zpisc(jpnpi,1) =  6.64e-7_wp
+      zpisc(jpppi,1) =  6.64e-7_wp
+      zpisc(jppfe,1) =  3.89e-11_wp
+      zpisc(jppch,1) =  1.17e-7_wp
+      zpisc(jpnph,1) =  6.64e-7_wp
+      zpisc(jppph,1) =  6.64e-7_wp
+      zpisc(jpndi,1) =  3.41e-7_wp
+      zpisc(jppdi,1) =  3.41e-7_wp
+      zpisc(jppon,1) =  4.84e-7_wp
+      zpisc(jppop,1) =  4.84e-7_wp
+      zpisc(jpdon,1) =  1.06e-5_wp
+      zpisc(jpdop,1) =  1.06e-5_wp
+      zpisc(jpgon,1) =  1.05e-8_wp
+      zpisc(jpgop,1) =  1.05e-8_wp
+#endif
+
  
       DO jn = jp_pcs0, jp_pcs1

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90

@@ -10 +10 @@
    !!             2.0  !  2007-12  (C. Ethe, G. Madec) from trcini.pisces.h90
    !!             3.5  !  2012-05  (C. Ethe) Merge PISCES-LOBSTER
-   !!----------------------------------------------------------------------
-#if defined key_pisces || defined key_pisces_reduced
-   !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
+   !!----------------------------------------------------------------------
+#if defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota
+   !!----------------------------------------------------------------------
+   !!   'key_pisces*'                                       PISCES bio-model
    !!----------------------------------------------------------------------
    !! trc_ini_pisces   : PISCES biochemical model initialisation
@@ -43 +44 @@
       !!----------------------------------------------------------------------
 
-      IF( lk_p4z ) THEN  ;   CALL p4z_ini   !  PISCES
-      ELSE               ;   CALL p2z_ini   !  LOBSTER
-      ENDIF
+      SELECT CASE ( nn_p4z )
+      !
+        CASE(1)          ;   CALL p2z_ini   !  LOBSTER
+        CASE(2)          ;   CALL p4z_ini   !  PISCES
+        CASE(3)          ;   CALL p5z_ini   !  PISCES QUOTA
+
+      END SELECT
 
    END SUBROUTINE trc_ini_pisces
+
+   SUBROUTINE p5z_ini
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE p5z_ini ***
+      !!
+      !! ** Purpose :   Initialisation of the PISCES biochemical model
+      !!                with variable stoichiometry
+      !!----------------------------------------------------------------------
+#if defined key_pisces_quota 
+      !
+      USE p5zsms          ! Main P4Z routine
+      USE p4zche          !  Chemical model
+      USE p5zsink         !  vertical flux of particulate matter due to sinking
+      USE p4zopt          !  optical model
+      USE p4zsbc          !  Boundary conditions
+      USE p4zfechem       !  Iron chemistry
+      USE p5zrem          !  Remineralisation of organic matter
+      USE p5zpoc          !  Remineralization of organic particles
+      USE p4zligand       !  Remineralization of organic ligands
+      USE p4zflx          !  Gas exchange
+      USE p5zlim          !  Co-limitations of differents nutrients
+      USE p5zprod         !  Growth rate of the 2 phyto groups
+      USE p5zmicro        !  Sources and sinks of microzooplankton
+      USE p5zmeso         !  Sources and sinks of mesozooplankton
+      USE p5zmort         !  Mortality terms for phytoplankton
+      USE p4zlys          !  Calcite saturation
+      USE p5zsed          !  Sedimentation & burial
+      !
+      REAL(wp), SAVE :: sco2   =  2.312e-3_wp
+      REAL(wp), SAVE :: alka0  =  2.423e-3_wp
+      REAL(wp), SAVE :: oxyg0  =  177.6e-6_wp 
+      REAL(wp), SAVE :: po4    =  2.174e-6_wp 
+      REAL(wp), SAVE :: bioma0 =  1.000e-8_wp  
+      REAL(wp), SAVE :: silic1 =  91.65e-6_wp  
+      REAL(wp), SAVE :: no3    =  31.04e-6_wp * 7.625_wp
+      !
+      INTEGER  ::  ierr
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' p5z_ini :   PISCES biochemical model initialisation'
+      IF(lwp) WRITE(numout,*) '             With variable stoichiometry'
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+                                                 ! Allocate PISCES arrays
+      ierr =         sms_pisces_alloc()          
+      ierr = ierr +  p4z_che_alloc()
+      ierr = ierr +  p5z_sink_alloc()
+      ierr = ierr +  p4z_opt_alloc()
+      ierr = ierr +  p5z_lim_alloc()
+      ierr = ierr +  p5z_prod_alloc()
+      ierr = ierr +  p5z_rem_alloc()
+      ierr = ierr +  p4z_flx_alloc()
+      ierr = ierr +  p5z_sed_alloc()
+      !
+      IF( lk_mpp    )   CALL mpp_sum( ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'pisces_alloc: unable to allocate PISCES arrays' )
+      !
+      ryyss    = nyear_len(1) * rday    ! number of seconds per year
+      r1_ryyss = 1. / ryyss
+      !
+      CALL p5z_sms_init       !  Maint routine
+      !                                            ! Time-step
+      ! Set biological ratios
+      ! ---------------------
+      rno3    =  16._wp / 122._wp
+      po4r    =   1._wp / 122._wp
+      o2nit   =  32._wp / 122._wp
+      rdenit  = 105._wp /  16._wp
+      rdenita =   3._wp /  5._wp
+      o2ut    = 133._wp / 122._wp
+      no3rat3 = no3rat3 / rno3
+      po4rat3 = po4rat3 / po4r
+
+      ! Initialization of tracer concentration in case of  no restart 
+      !--------------------------------------------------------------
+      IF( .NOT. ln_rsttr ) THEN  
+         
+         trn(:,:,:,jpdic) = sco2
+         trn(:,:,:,jpdoc) = bioma0
+         trn(:,:,:,jpdon) = bioma0
+         trn(:,:,:,jpdop) = bioma0
+         trn(:,:,:,jptal) = alka0
+         trn(:,:,:,jpoxy) = oxyg0
+         trn(:,:,:,jpcal) = bioma0
+         trn(:,:,:,jppo4) = po4 / po4r
+         trn(:,:,:,jppoc) = bioma0
+         trn(:,:,:,jppon) = bioma0
+         trn(:,:,:,jppop) = bioma0
+#  if ! defined key_kriest
+         trn(:,:,:,jpgoc) = bioma0
+         trn(:,:,:,jpgon) = bioma0
+         trn(:,:,:,jpgop) = bioma0
+         trn(:,:,:,jpbfe) = bioma0 * 5.e-6
+#  else
+         trn(:,:,:,jpnum) = bioma0 / ( 6. * xkr_massp )
+#  endif
+         trn(:,:,:,jpsil) = silic1
+         trn(:,:,:,jpdsi) = bioma0 * 0.15
+         trn(:,:,:,jpgsi) = bioma0 * 5.e-6
+         trn(:,:,:,jpphy) = bioma0
+         trn(:,:,:,jpnph) = bioma0
+         trn(:,:,:,jppph) = bioma0
+         trn(:,:,:,jppic) = bioma0
+         trn(:,:,:,jpnpi) = bioma0
+         trn(:,:,:,jpppi) = bioma0
+         trn(:,:,:,jpdia) = bioma0
+         trn(:,:,:,jpndi) = bioma0
+         trn(:,:,:,jppdi) = bioma0
+         trn(:,:,:,jpzoo) = bioma0
+         trn(:,:,:,jpmes) = bioma0
+         trn(:,:,:,jpfer) = 0.6E-9
+         trn(:,:,:,jpsfe) = bioma0 * 5.e-6
+         trn(:,:,:,jppfe) = bioma0 * 5.e-6
+         trn(:,:,:,jpdfe) = bioma0 * 5.e-6
+         trn(:,:,:,jpnfe) = bioma0 * 5.e-6
+         trn(:,:,:,jpnch) = bioma0 * 12. / 55.
+         trn(:,:,:,jppch) = bioma0 * 12. / 55.
+         trn(:,:,:,jpdch) = bioma0 * 12. / 55.
+         trn(:,:,:,jpno3) = no3
+         trn(:,:,:,jpnh4) = bioma0
+#if defined key_ligand
+         trn(:,:,:,jplgw) = 0.6E-9
+         trn(:,:,:,jpfep) = 0. * 5.e-6
+#endif
+
+         ! initialize the half saturation constant for silicate
+         ! ----------------------------------------------------
+         xksi(:,:)    = 2.e-6
+         xksimax(:,:) = xksi(:,:)
+      END IF
+
+      CALL p5z_sink_init      !  vertical flux of particulate organic matter
+      CALL p4z_opt_init       !  Optic: PAR in the water column
+      CALL p5z_lim_init       !  co-limitations by the various nutrients
+      CALL p5z_prod_init      !  phytoplankton growth rate over the global ocean.
+      CALL p4z_sbc_init       !  boundary conditions
+      CALL p4z_fechem_init    !  Iron chemistry
+      CALL p5z_rem_init       !  remineralisation
+      CALL p5z_poc_init       !  remineralisation of organic particles
+#if defined key_ligand
+      CALL p4z_ligand_init    !  remineralisation of organic ligands
+#endif
+      CALL p5z_mort_init      !  phytoplankton mortality 
+      CALL p5z_micro_init     !  microzooplankton
+      CALL p5z_meso_init      !  mesozooplankton
+      CALL p4z_lys_init       !  calcite saturation
+      CALL p4z_flx_init       !  gas exchange 
+
+      ndayflxtr = 0
+
+      IF(lwp) WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*) 'Initialization of PISCES tracers done'
+      IF(lwp) WRITE(numout,*) 
+#endif
+      !
+   END SUBROUTINE p5z_ini
 
    SUBROUTINE p4z_ini
@@ -64 +226 @@
       USE p4zfechem       !  Iron chemistry
       USE p4zrem          !  Remineralisation of organic matter
+      USE p4zpoc          !  Remineralisation of organic particles
+      USE p4zligand       !  Remineralization of organic ligands
       USE p4zflx          !  Gas exchange
       USE p4zlim          !  Co-limitations of differents nutrients
@@ -74 +238 @@
       !
       REAL(wp), SAVE :: sco2   =  2.312e-3_wp
-      REAL(wp), SAVE :: alka0  =  2.426e-3_wp
-      REAL(wp), SAVE :: oxyg0  =  177.6e-6_wp 
-      REAL(wp), SAVE :: po4    =  2.165e-6_wp 
-      REAL(wp), SAVE :: bioma0 =  1.000e-8_wp  
-      REAL(wp), SAVE :: silic1 =  91.51e-6_wp  
-      REAL(wp), SAVE :: no3    =  30.9e-6_wp * 7.625_wp
-      !
-      INTEGER  ::  ji, jj, jk, ierr
-      REAL(wp) ::  zcaralk, zbicarb, zco3
-      REAL(wp) ::  ztmas, ztmas1
+      REAL(wp), SAVE :: alka0  =  2.423e-3_wp
+      REAL(wp), SAVE :: oxyg0  =  177.6e-6_wp
+      REAL(wp), SAVE :: po4    =  2.174e-6_wp
+      REAL(wp), SAVE :: bioma0 =  1.000e-8_wp
+      REAL(wp), SAVE :: silic1 =  91.65e-6_wp
+      REAL(wp), SAVE :: no3    =  31.04e-6_wp * 7.625_wp
+      !
+      INTEGER  ::  ierr
       !!----------------------------------------------------------------------
 
@@ -91 +253 @@
 
                                                  ! Allocate PISCES arrays
-      ierr =         sms_pisces_alloc()          
+      ierr =         sms_pisces_alloc()
       ierr = ierr +  p4z_che_alloc()
       ierr = ierr +  p4z_sink_alloc()
       ierr = ierr +  p4z_opt_alloc()
+      ierr = ierr +  p4z_lim_alloc()
       ierr = ierr +  p4z_prod_alloc()
       ierr = ierr +  p4z_rem_alloc()
@@ -109 +272 @@
       CALL p4z_sms_init       !  Maint routine
       !                                            ! Time-step
-
       ! Set biological ratios
       ! ---------------------
@@ -121 +283 @@
       ! Initialization of tracer concentration in case of  no restart 
       !--------------------------------------------------------------
-      IF( .NOT. ln_rsttr ) THEN  
-         
+      IF( .NOT. ln_rsttr ) THEN
+
          trn(:,:,:,jpdic) = sco2
          trn(:,:,:,jpdoc) = bioma0
@@ -151 +313 @@
          trn(:,:,:,jpno3) = no3
          trn(:,:,:,jpnh4) = bioma0
+#if defined key_ligand
+         trn(:,:,:,jplgw) = 0.6E-9
+         trn(:,:,:,jpfep) = 0. * 5.e-6
+#endif
 
          ! initialize the half saturation constant for silicate
@@ -157 +323 @@
          xksimax(:,:) = xksi(:,:)
       END IF
-
 
       CALL p4z_sink_init      !  vertical flux of particulate organic matter
@@ -166 +331 @@
       CALL p4z_fechem_init    !  Iron chemistry
       CALL p4z_rem_init       !  remineralisation
+      CALL p4z_poc_init       !  remineralization of organic particles
+#if defined key_ligand
+      CALL p4z_ligand_init    !  remineralisation of organic ligands
+#endif
       CALL p4z_mort_init      !  phytoplankton mortality 
       CALL p4z_micro_init     !  microzooplankton
@@ -174 +343 @@
       ndayflxtr = 0
 
-      IF(lwp) WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'Initialization of PISCES tracers done'
-      IF(lwp) WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*)
 #endif
       !

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcnam_pisces.F90

@@ -8 +8 @@
    !!             1.0  !  2003-08 (C. Ethe)  module F90
    !!             2.0  !  2007-12  (C. Ethe, G. Madec) from trcnam.pisces.h90
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
-#if defined key_pisces || defined key_pisces_reduced
+#if defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota
    !!----------------------------------------------------------------------
-   !!   'key_pisces'   :                                   PISCES bio-model
+   !!   'key_pisces*'  :                                   PISCES bio-model
    !!----------------------------------------------------------------------
    !! trc_nam_pisces       : PISCES model namelist read
@@ -66 +67 @@
 #if defined key_pisces
       IF(lwp) WRITE(numout,*) ' trc_nam_pisces : read PISCES namelist'
+#elif defined key_pisces_quota
+      IF(lwp) WRITE(numout,*) ' trc_nam_pisces : read PISCES-QUOTA namelist'
 #else
       IF(lwp) WRITE(numout,*) ' trc_nam_pisces : read LOBSTER namelist'
@@ -123 +126 @@
 #if defined key_pisces_reduced
 
-      IF( ( .NOT.lk_iomput .AND. ln_diabio ) .OR. lk_trdmxl_trc ) THEN
+      IF( ( .NOT.lk_iomput .AND. ln_diabio ) .OR. lk_trdmld_trc ) THEN
          !
          ! Namelist nampisdbi

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcsms_pisces.F90

@@ -6 +6 @@
    !! History :   1.0  !  2004-03 (O. Aumont) Original code
    !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
-#if defined key_pisces || defined key_pisces_reduced
+#if defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota
    !!----------------------------------------------------------------------
-   !!   'key_pisces'                                       PISCES bio-model
+   !!   'key_pisces*'                                       PISCES bio-model
    !!----------------------------------------------------------------------
    !!   trcsms_pisces        :  Time loop of passive tracers sms
    !!----------------------------------------------------------------------
    USE par_pisces
+   USE p5zsms
    USE p4zsms
    USE p2zsms
@@ -48 +50 @@
       !!---------------------------------------------------------------------
       !
-      IF( lk_p4z ) THEN  ;   CALL p4z_sms( kt )   !  PISCES
-      ELSE               ;   CALL p2z_sms( kt )   !  LOBSTER
-      ENDIF
+      SELECT CASE ( nn_p4z )
+
+      CASE(1)       ;        CALL p2z_sms( kt )   !  LOBSTER
+      CASE(2)       ;        CALL p4z_sms( kt )   !  PISCES
+      CASE(3)       ;        CALL p5z_sms( kt )   !  PISCES with variable stoichiometry
+
+      END SELECT
 
       !

branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90

@@ -5 +5 @@
    !!======================================================================
    !! History :   1.0  !  2009-05 (C. Ethe)  Original code
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
-#if defined key_top && defined key_iomput && ( defined key_pisces || defined key_pisces_reduced )
+#if defined key_top && defined key_iomput && ( defined key_pisces || defined key_pisces_reduced || defined key_pisces_quota)
    !!----------------------------------------------------------------------
-   !!   'key_pisces or key_pisces_reduced'                     PISCES model
+   !!   'key_pisces*'                     PISCES model
    !!----------------------------------------------------------------------
    !! trc_wri_pisces   :  outputs of concentration fields
@@ -30 +31 @@
       !! ** Purpose :   output passive tracers fields 
       !!---------------------------------------------------------------------
-      CHARACTER (len=20)           :: cltra
-      REAL(wp)                     :: zfact
-      INTEGER                      :: ji, jj, jk, jn
-      REAL(wp), DIMENSION(jpi,jpj) :: zdic, zo2min, zdepo2min
+      CHARACTER (len=20)   :: cltra
+      REAL(wp)             :: zrfact
+      INTEGER              :: jn
       !!---------------------------------------------------------------------
  
@@ -41 +41 @@
       DO jn = jp_pcs0, jp_pcs1
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         CALL iom_put( cltra, trn(:,:,:,jn) )
+         IF( lk_vvl ) THEN
+            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) )
+         ELSE
+            CALL iom_put( cltra, trn(:,:,:,jn) )
+         ENDIF
+         CALL iom_put( cltra, trn(:,:,:,jn) * zrfact )
       END DO
 #else
       DO jn = jp_pcs0, jp_pcs1
-         zfact = 1.0e+6 
-         IF( jn == jpno3 .OR. jn == jpnh4 ) zfact = rno3 * 1.0e+6 
-         IF( jn == jppo4  )                 zfact = po4r * 1.0e+6
+         zrfact = 1.0e+6 
+         IF( jn == jpno3 .OR. jn == jpnh4 ) zrfact = rno3 * 1.0e+6 
+         IF( jn == jppo4  )                 zrfact = po4r * 1.0e+6
          cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer
-         IF( iom_use( cltra ) )  CALL iom_put( cltra, trn(:,:,:,jn) * zfact )
+         IF( lk_vvl ) THEN
+            CALL iom_put( cltra, trn(:,:,:,jn) * fse3t_n(:,:,:) * zrfact )
+         ELSE
+            CALL iom_put( cltra, trn(:,:,:,jn) * zrfact )
+         ENDIF
       END DO
-
-      IF( iom_use( "INTDIC" ) ) THEN                     !   DIC content in kg/m2
-         zdic(:,:) = 0.
-         DO jk = 1, jpkm1
-            zdic(:,:) = zdic(:,:) + trn(:,:,jk,jpdic) * fse3t(:,:,jk) * tmask(:,:,jk) * 12.
-         ENDDO
-         CALL iom_put( 'INTDIC', zdic )     
-      ENDIF
-      !
-      IF( iom_use( "O2MIN" ) .OR. iom_use ( "ZO2MIN" ) ) THEN  ! Oxygen minimum concentration and depth 
-         zo2min   (:,:) = trn(:,:,1,jpoxy) * tmask(:,:,1)
-         zdepo2min(:,:) = fsdepw(:,:,1)    * tmask(:,:,1)
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( tmask(ji,jj,jk) == 1 ) then
-                     IF( trn(ji,jj,jk,jpoxy) < zo2min(ji,jj) ) then
-                        zo2min   (ji,jj) = trn(ji,jj,jk,jpoxy)
-                        zdepo2min(ji,jj) = fsdepw(ji,jj,jk)
-                     ENDIF
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-         CALL iom_put('O2MIN' , zo2min     )                              ! oxygen minimum concentration
-         CALL iom_put('ZO2MIN', zdepo2min  )                              ! depth of oxygen minimum concentration
-          !
-      ENDIF
 #endif
       !