New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9987 for branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

Ignore:
Timestamp:
2018-07-23T11:33:03+02:00 (6 years ago)
Author:
emmafiedler
Message:

Merge with GO6 FOAMv14 package branch r9288

Location:
branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z
Files:
16 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90

    r7959 r9987  
    109109 
    110110   !!====================================================================== 
    111 END MODULE  p4zbio 
    112  
     111END MODULE p4zbio 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90

    r7960 r9987  
    3232   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fekeq    ! chemistry of Fe 
    3333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemc    ! Solubilities of O2 and CO2 
    34    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2    ! Solubilities of O2 and CO2 
     34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2   ! Solubilities of O2 and CO2 
     35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature 
    3536 
    3637   REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm 
     
    3940   REAL(wp) ::   o2atm  = 1. / ( 1000. * 0.20946 )   
    4041 
    41    REAL(wp) ::   akcc1  = -171.9065       ! coeff. for apparent solubility equilibrium 
    42    REAL(wp) ::   akcc2  =   -0.077993     ! Millero et al. 1995 from Mucci 1983 
    43    REAL(wp) ::   akcc3  = 2839.319         
    44    REAL(wp) ::   akcc4  =   71.595         
    45    REAL(wp) ::   akcc5  =   -0.77712       
    46    REAL(wp) ::   akcc6  =    0.00284263    
    47    REAL(wp) ::   akcc7  =  178.34         
    48    REAL(wp) ::   akcc8  =   -0.07711      
    49    REAL(wp) ::   akcc9  =    0.0041249    
    50  
    51    REAL(wp) ::   rgas   = 83.143         ! universal gas constants 
     42   REAL(wp) ::   rgas   = 83.14472       ! universal gas constants 
    5243   REAL(wp) ::   oxyco  = 1. / 22.4144   ! converts from liters of an ideal gas to moles 
    5344 
     
    5546   REAL(wp) ::   bor2   = 1. / 10.82 
    5647 
    57    REAL(wp) ::   ca0    = -162.8301      ! WEISS & PRICE 1980, units mol/(kg atm) 
    58    REAL(wp) ::   ca1    =  218.2968 
    59    REAL(wp) ::   ca2    =   90.9241 
    60    REAL(wp) ::   ca3    =   -1.47696 
    61    REAL(wp) ::   ca4    =    0.025695 
    62    REAL(wp) ::   ca5    =   -0.025225 
    63    REAL(wp) ::   ca6    =    0.0049867 
    64  
    65    REAL(wp) ::   c10    = -3670.7        ! Coeff. for 1. dissoc. of carbonic acid (Edmond and Gieskes, 1970)    
    66    REAL(wp) ::   c11    =    62.008      
    67    REAL(wp) ::   c12    =    -9.7944     
    68    REAL(wp) ::   c13    =     0.0118      
    69    REAL(wp) ::   c14    =    -0.000116 
    70  
    71    REAL(wp) ::   c20    = -1394.7       ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)    
    72    REAL(wp) ::   c21    =    -4.777    
    73    REAL(wp) ::   c22    =     0.0184    
    74    REAL(wp) ::   c23    =    -0.000118 
    75  
    7648   REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate 
    7749   REAL(wp) ::   st2    =  1./96.062    !  (Morris & Riley 1966) 
    78    REAL(wp) ::   ks0    =    141.328  
    79    REAL(wp) ::   ks1    =  -4276.1   
    80    REAL(wp) ::   ks2    =    -23.093 
    81    REAL(wp) ::   ks3    = -13856.   
    82    REAL(wp) ::   ks4    =   324.57  
    83    REAL(wp) ::   ks5    =   -47.986 
    84    REAL(wp) ::   ks6    =  35474.  
    85    REAL(wp) ::   ks7    =   -771.54 
    86    REAL(wp) ::   ks8    =    114.723 
    87    REAL(wp) ::   ks9    =  -2698.   
    88    REAL(wp) ::   ks10   =   1776.  
    89    REAL(wp) ::   ks11   =      1. 
    90    REAL(wp) ::   ks12   =     -0.001005  
    9150 
    9251   REAL(wp) ::   ft1    =    0.000067   ! constants for calculate concentrations for fluorides 
    9352   REAL(wp) ::   ft2    = 1./18.9984    ! (Dickson & Riley 1979 ) 
    94    REAL(wp) ::   kf0    =  -12.641     
    95    REAL(wp) ::   kf1    = 1590.2     
    96    REAL(wp) ::   kf2    =    1.525     
    97    REAL(wp) ::   kf3    =    1.0      
    98    REAL(wp) ::   kf4    =   -0.001005 
    99  
    100    REAL(wp) ::   cb0    = -8966.90      ! Coeff. for 1. dissoc. of boric acid  
    101    REAL(wp) ::   cb1    = -2890.53      ! (Dickson and Goyet, 1994) 
    102    REAL(wp) ::   cb2    =   -77.942 
    103    REAL(wp) ::   cb3    =     1.728 
    104    REAL(wp) ::   cb4    =    -0.0996 
    105    REAL(wp) ::   cb5    =   148.0248 
    106    REAL(wp) ::   cb6    =   137.1942 
    107    REAL(wp) ::   cb7    =     1.62142 
    108    REAL(wp) ::   cb8    =   -24.4344 
    109    REAL(wp) ::   cb9    =   -25.085 
    110    REAL(wp) ::   cb10   =    -0.2474  
    111    REAL(wp) ::   cb11   =     0.053105 
    112  
    113    REAL(wp) ::   cw0    = -13847.26     ! Coeff. for dissoc. of water (Dickson and Riley, 1979 ) 
    114    REAL(wp) ::   cw1    =    148.9652   
    115    REAL(wp) ::   cw2    =    -23.6521 
    116    REAL(wp) ::   cw3    =    118.67  
    117    REAL(wp) ::   cw4    =     -5.977  
    118    REAL(wp) ::   cw5    =      1.0495   
    119    REAL(wp) ::   cw6    =     -0.01615 
    12053 
    12154   !                                    ! volumetric solubility constants for o2 in ml/L   
     
    185118      REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 
    186119      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2 
    187       REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1 
    188       REAL(wp) ::   zis  , zis2 , zsal15, zisqrt 
     120      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat 
     121      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1  , za2 
    189122      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw 
    190123      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1 
     
    193126      IF( nn_timing == 1 )  CALL timing_start('p4z_che') 
    194127      ! 
     128      ! Computations of chemical constants require in situ temperature 
     129      ! Here a quite simple formulation is used to convert  
     130      ! potential temperature to in situ temperature. The errors is less than  
     131      ! 0.04°C relative to an exact computation 
     132      ! --------------------------------------------------------------------- 
     133      DO jk = 1, jpk 
     134         DO jj = 1, jpj 
     135            DO ji = 1, jpi 
     136               zpres = fsdept(ji,jj,jk) / 1000. 
     137               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (tsn(ji,jj,jk,jp_sal) - 35.0) ) 
     138               za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 ) 
     139               tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2 
     140            END DO 
     141         END DO 
     142      END DO 
     143      ! 
    195144      ! CHEMICAL CONSTANTS - SURFACE LAYER 
    196145      ! ---------------------------------- 
     
    200149         DO ji = 1, jpi 
    201150            !                             ! SET ABSOLUTE TEMPERATURE 
    202             ztkel = tsn(ji,jj,1,jp_tem) + 273.16 
     151            ztkel = tempis(ji,jj,1) + 273.15 
    203152            zt    = ztkel * 0.01 
    204153            zt2   = zt * zt 
     
    208157            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 
    209158            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
    210             zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 ) 
    211             !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 in ml/L (EQ. 8, GARCIA AND GORDON, 1992) 
    212             ztgg  = LOG( ( 298.15 - tsn(ji,jj,1,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
    213             ztgg2 = ztgg  * ztgg 
    214             ztgg3 = ztgg2 * ztgg 
    215             ztgg4 = ztgg3 * ztgg 
    216             ztgg5 = ztgg4 * ztgg 
    217             zoxy  = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5   & 
    218                    + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) +  ox10 * zsal2 
    219  
     159            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    & 
     160            &       + 0.0047036e-4*ztkel**2) 
    220161            !                             ! SET SOLUBILITIES OF O2 AND CO2  
    221             chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm) 
    222             chemc(ji,jj,2) = ( EXP( zoxy  ) * o2atm ) * oxyco              ! mol/(L atm) 
     162            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm) 
     163            chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 
     164            chemc(ji,jj,3) = 57.7 - 0.118*ztkel 
    223165            ! 
    224166         END DO 
     
    233175!CDIR NOVERRCHK 
    234176            DO ji = 1, jpi 
    235               ztkel = tsn(ji,jj,jk,jp_tem) + 273.16 
     177              ztkel = tempis(ji,jj,jk) + 273.15 
    236178              zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35. 
    237179              zsal2 = zsal * zsal 
    238               ztgg  = LOG( ( 298.15 - tsn(ji,jj,jk,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
     180              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
    239181              ztgg2 = ztgg  * ztgg 
    240182              ztgg3 = ztgg2 * ztgg 
     
    259201            DO ji = 1, jpi 
    260202 
    261                ! SET PRESSION 
    262                zpres   = 1.025e-1 * fsdept(ji,jj,jk) 
     203               ! SET PRESSION ACCORDING TO SAUNDER (1980) 
     204               zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) ) 
     205               zc1 = 5.92E-3 + zplat**2 * 5.25E-3 
     206               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*fsdept(ji,jj,jk)))) / 4.42E-6 
     207               zpres = zpres / 10.0 
    263208 
    264209               ! SET ABSOLUTE TEMPERATURE 
    265                ztkel   = tsn(ji,jj,jk,jp_tem) + 273.16 
     210               ztkel   = tempis(ji,jj,jk) + 273.15 
    266211               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35. 
    267212               zsqrt  = SQRT( zsal ) 
     
    272217               zis2   = zis * zis 
    273218               zisqrt = SQRT( zis ) 
    274                ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20. 
     219               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20. 
    275220 
    276221               ! CHLORINITY (WOOSTER ET AL., 1969) 
     
    284229 
    285230               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 
    286                zcks    = EXP(  ks1 * ztr + ks0 + ks2 * zlogt                           & 
    287                   &                     + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt   & 
    288                   &                     + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis      & 
    289                   &                     + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal )  ) 
     231               zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         & 
     232               &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt & 
     233               &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    & 
     234               &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         & 
     235               &         + LOG(1.0 - 0.001005 * zsal)) 
     236               ! 
     237               aphscale(ji,jj,jk) = ( 1. + zst / zcks ) 
    290238 
    291239               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 
    292                zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  ) 
     240               zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   & 
     241               &         + LOG(1.0d0 - 0.001005d0*zsal)            & 
     242               &         + LOG(1.0d0 + zst/zcks)) 
    293243 
    294244               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 
    295                zckb    = ( cb0 + cb1 * zsqrt + cb2  * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr   & 
    296                   &    + ( cb5 + cb6 * zsqrt + cb7  * zsal )                                            & 
    297                   &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             & 
    298                   &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  ) 
    299  
    300                zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal 
    301                zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2 
     245               zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        & 
     246               &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         & 
     247               &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   & 
     248               &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      &  
     249               &      * zlogt + 0.053105*zsqrt*ztkel 
     250 
     251 
     252               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO  
     253               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 
     254               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  & 
     255                  - 0.011555*zsal + 0.0001152*zsal*zsal) 
     256               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      & 
     257                  - 0.01781*zsal + 0.0001122*zsal*zsal) 
    302258 
    303259               ! PKW (H2O) (DICKSON AND RILEY, 1979) 
    304                zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal 
    305  
     260               zckw = -13847.26*ztr + 148.9652 - 23.6521 * zlogt    &  
     261               &     + (118.67*ztr - 5.977 + 1.0495 * zlogt)        & 
     262               &     * zsqrt - 0.01615 * zsal 
    306263 
    307264               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
    308265               !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 
    309                zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   & 
    310                   &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15 
     266               zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   & 
     267                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  & 
     268                  &      - 0.07711*zsal + 0.0041249*zsal15 
    311269 
    312270               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 
     
    378336      !!                     ***  ROUTINE p4z_che_alloc  *** 
    379337      !!---------------------------------------------------------------------- 
    380       ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,2), chemo2(jpi,jpj,jpk), STAT=p4z_che_alloc ) 
     338      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk),   & 
     339      &         tempis(jpi,jpj,jpk), STAT=p4z_che_alloc ) 
    381340      ! 
    382341      IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') 
     
    396355 
    397356   !!====================================================================== 
    398 END MODULE  p4zche 
     357END MODULE p4zche 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90

    r7960 r9987  
    8484      ! 
    8585      INTEGER  ::   ji, jj, jm, iind, iindm1 
    86       REAL(wp) ::   ztc, ztc2, ztc3, zws, zkgwan 
     86      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan 
    8787      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact 
     88      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff 
    8889      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 
    8990      REAL(wp) ::   zyr_dec, zdco2dt 
    9091      CHARACTER (len=25) :: charout 
    91       REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d  
     92      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d, zpco2atm  
    9293      !!--------------------------------------------------------------------- 
    9394      ! 
    9495      IF( nn_timing == 1 )  CALL timing_start('p4z_flx') 
    9596      ! 
    96       CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx ) 
     97      CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 
    9798      ! 
    9899 
     
    135136 
    136137               ! CALCULATE [ALK]([CO3--], [HCO3-]) 
    137                zalk  = zalka - (  akw3(ji,jj,1) / zph - zph + zbot / ( 1.+ zph / akb3(ji,jj,1) )  ) 
     138               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph / aphscale(ji,jj,1)    & 
     139               &       + zbot / ( 1.+ zph / akb3(ji,jj,1) )  ) 
    138140 
    139141               ! CALCULATE [H+] AND [H2CO3] 
     
    162164            ztc2 = ztc * ztc 
    163165            ztc3 = ztc * ztc2  
     166            ztc4 = ztc2 * ztc2  
    164167            ! Compute the schmidt Number both O2 and CO2 
    165             zsch_co2 = 2073.1 - 125.62 * ztc + 3.6276 * ztc2 - 0.043126 * ztc3 
    166             zsch_o2  = 1953.4 - 128.0  * ztc + 3.9918 * ztc2 - 0.050091 * ztc3 
     168            zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4 
     169            zsch_o2  = 1920.4 - 135.6  * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4 
    167170            !  wind speed  
    168171            zws  = wndm(ji,jj) * wndm(ji,jj) 
    169172            ! Compute the piston velocity for O2 and CO2 
    170             zkgwan = 0.3 * zws  + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946  * ztc2 ) 
     173            zkgwan = 0.251 * zws 
    171174            zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 
    172175# if defined key_degrad 
     
    181184      DO jj = 1, jpj 
    182185         DO ji = 1, jpi 
     186            ztkel  = tsn(ji,jj,1,jp_tem) + 273.15 
     187            zsal   = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35. 
     188            zvapsw = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal) 
     189            zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) 
     190            zxc2 = (1.0 - zpco2atm(ji,jj) * 1E-6 )**2 
     191            zfugcoeff = EXP(patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   & 
     192            &           / (82.05736 * ztkel)) 
     193            zfco2 = zpco2atm(ji,jj) * zfugcoeff 
     194 
    183195            ! Compute CO2 flux for the sea and air 
    184             zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)   ! (mol/L) * (m/s) 
    185             zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ? 
     196            zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s) 
     197            zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ? 
    186198            oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 
    187199            ! compute the trend 
    188             tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1) 
     200            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1) * tmask(ji,jj,1) 
    189201 
    190202            ! Compute O2 flux  
    191             zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
    192             zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj) 
    193             zoflx(ji,jj) = zfld16 - zflu16 
     203            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
     204            zflu16 = trb(ji,jj,1,jpoxy) * zkgo2(ji,jj) 
     205            zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1) 
    194206            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / fse3t(ji,jj,1) 
    195207         END DO 
     
    222234         ENDIF 
    223235         IF( iom_use( "Dpco2" ) ) THEN 
    224            zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
     236           zw2d(:,:) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
    225237           CALL iom_put( "Dpco2" ,  zw2d ) 
    226238         ENDIF 
    227239         IF( iom_use( "Dpo2" ) )  THEN 
    228            zw2d(:,:) = ( atcox * patm(:,:) - trb(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1) 
     240           zw2d(:,:) = ( atcox * patm(:,:) - atcox * trn(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) 
    229241           CALL iom_put( "Dpo2"  , zw2d ) 
    230242         ENDIF 
     
    238250            trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1)  
    239251            trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1)  
    240             trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)  
    241          ENDIF 
    242       ENDIF 
    243       ! 
    244       CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx ) 
     252            trc2d(:,:,jp_pcs0_2d + 3) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
     253         ENDIF 
     254      ENDIF 
     255      ! 
     256      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 
    245257      ! 
    246258      IF( nn_timing == 1 )  CALL timing_stop('p4z_flx') 
     
    400412 
    401413   !!====================================================================== 
    402 END MODULE  p4zflx 
     414END MODULE p4zflx 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90

    r7960 r9987  
    8181 
    8282   !!====================================================================== 
    83 END MODULE  p4zint 
     83END MODULE p4zint 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90

    r7959 r9987  
    4444   REAL(wp), PUBLIC ::  xkdoc       !:  2nd half-sat. of DOC remineralization   
    4545   REAL(wp), PUBLIC ::  concbfe     !:  Fe half saturation for bacteria  
     46   REAL(wp), PUBLIC ::  oxymin      !:  half saturation constant for anoxia 
    4647   REAL(wp), PUBLIC ::  qnfelim     !:  optimal Fe quota for nanophyto 
    4748   REAL(wp), PUBLIC ::  qdfelim     !:  optimal Fe quota for diatoms 
     
    121122               zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) 
    122123               zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concbnh4 ) 
    123                zlim3    = trb(ji,jj,jk,jpfer) / ( concbfe + trb(ji,jj,jk,jpfer) ) 
     124               zlim3    = biron(ji,jj,jk)     / ( concbfe + biron(ji,jj,jk) ) 
    124125               zlim4    = trb(ji,jj,jk,jpdoc) / ( xkdoc   + trb(ji,jj,jk,jpdoc) ) 
    125126               xlimbacl(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 ) 
     
    187188      END DO 
    188189      ! 
     190      DO jk = 1, jpkm1 
     191         DO jj = 1, jpj 
     192            DO ji = 1, jpi 
     193               ! denitrification factor computed from O2 levels 
     194               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    & 
     195                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  ) 
     196               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) ) 
     197            END DO 
     198         END DO 
     199      END DO 
    189200      ! 
    190201      IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics 
     
    216227      NAMELIST/nampislim/ concnno3, concdno3, concnnh4, concdnh4, concnfer, concdfer, concbfe,   & 
    217228         &                concbno3, concbnh4, xsizedia, xsizephy, xsizern, xsizerd,          &  
    218          &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r 
     229         &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r, oxymin 
    219230      INTEGER :: ios                 ! Local integer output status for namelist read 
    220231 
     
    249260         WRITE(numout,*) '    Minimum size criteria for nanophyto      xsizephy  = ', xsizephy 
    250261         WRITE(numout,*) '    Fe half saturation for bacteria          concbfe   = ', concbfe 
     262         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin   =' , oxymin 
    251263         WRITE(numout,*) '    optimal Fe quota for nano.               qnfelim   = ', qnfelim 
    252264         WRITE(numout,*) '    Optimal Fe quota for diatoms             qdfelim   = ', qdfelim 
    253265      ENDIF 
    254  
     266      ! 
     267      nitrfac (:,:,:) = 0._wp 
     268      ! 
    255269   END SUBROUTINE p4z_lim_init 
    256270 
     
    265279 
    266280   !!====================================================================== 
    267 END MODULE  p4zlim 
     281END MODULE p4zlim 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

    r7959 r9987  
    6565      REAL(wp) ::   zomegaca, zexcess, zexcess0 
    6666      CHARACTER (len=25) :: charout 
    67       REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss    
     67      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss    
    6868      !!--------------------------------------------------------------------- 
    6969      ! 
    7070      IF( nn_timing == 1 )  CALL timing_start('p4z_lys') 
    7171      ! 
    72       CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss ) 
     72      CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss ) 
    7373      ! 
    7474      zco3    (:,:,:) = 0. 
     
    9191                  zalka = trb(ji,jj,jk,jptal) / zfact 
    9292                  ! CALCULATE [ALK]([CO3--], [HCO3-]) 
    93                   zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 
     93                  zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph / ( aphscale(ji,jj,jk) + rtrn )  & 
     94                  &       + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 
    9495                  ! CALCULATE [H+] and [CO3--] 
    9596                  zaldi = zdic - zalk 
     
    119120               zcalcon  = calcon * ( tsn(ji,jj,jk,jp_sal) / 35._wp ) 
    120121               zfact    = rhop(ji,jj,jk) / 1000._wp 
    121                zomegaca = ( zcalcon * zco3(ji,jj,jk) * zfact ) / aksp(ji,jj,jk)  
     122               zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn ) 
     123               zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn ) 
    122124 
    123125               ! SET DEGREE OF UNDER-/SUPERSATURATION 
     
    148150      IF( lk_iomput .AND. knt == nrdttrc ) THEN 
    149151         IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )          * tmask(:,:,:) ) 
    150          IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:) * 1.e+3               * tmask(:,:,:) ) 
    151          IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", aksp(:,:,:) * 1.e+3 / calcon      * tmask(:,:,:) ) 
    152          IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) ) 
     152         IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:)    * 1.e+3            * tmask(:,:,:) ) 
     153         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", zco3sat(:,:,:) * 1.e+3            * tmask(:,:,:) ) 
     154         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) ) 
    153155      ELSE 
    154          trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:) 
    155          trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:) 
    156          trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:) 
     156         IF( ln_diatrc ) THEN 
     157            trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:) 
     158            trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:) 
     159            trc3d(:,:,:,jp_pcs0_3d + 2) = zco3sat(:,:,:)           * tmask(:,:,:) 
     160         ENDIF 
    157161      ENDIF 
    158162      ! 
     
    163167      ENDIF 
    164168      ! 
    165       CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss ) 
     169      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss ) 
    166170      ! 
    167171      IF( nn_timing == 1 )  CALL timing_stop('p4z_lys') 
     
    223227#endif  
    224228   !!====================================================================== 
    225 END MODULE  p4zlys 
     229END MODULE p4zlys 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90

    r7959 r9987  
    340340 
    341341   !!====================================================================== 
    342 END MODULE  p4zmeso 
     342END MODULE p4zmeso 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90

    r7959 r9987  
    273273 
    274274   !!====================================================================== 
    275 END MODULE  p4zmicro 
     275END MODULE p4zmicro 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90

    r7960 r9987  
    277277 
    278278   !!====================================================================== 
    279 END MODULE  p4zmort 
     279END MODULE p4zmort 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90

    r7959 r9987  
    7676      REAL(wp) ::   zchl 
    7777      REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep 
    78       REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100 
     78      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 
     79      REAL(wp), POINTER, DIMENSION(:,:  ) :: zqsr100, zqsr_corr 
    7980      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3 
    8081      !!--------------------------------------------------------------------- 
     
    8384      ! 
    8485      ! Allocate temporary workspace 
    85       CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     86      CALL wrk_alloc( jpi, jpj,      zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     87      CALL wrk_alloc( jpi, jpj,      zqsr100, zqsr_corr ) 
    8688      CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 ) 
    8789 
     
    112114      !                                        !  -------------------------------------- 
    113115      IF( l_trcdm2dc ) THEN                     !  diurnal cycle 
    114          ! 1% of qsr to compute euphotic layer 
    115          zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr 
    116          ! 
    117          CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 )  
     116         ! 
     117         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
     118         ! 
     119         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )  
    118120         ! 
    119121         DO jk = 1, nksrp       
     
    123125         END DO 
    124126         ! 
    125          CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 )  
     127         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
     128         ! 
     129         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )  
    126130         ! 
    127131         DO jk = 1, nksrp       
     
    130134         ! 
    131135      ELSE 
    132          ! 1% of qsr to compute euphotic layer 
    133          zqsr100(:,:) = 0.01 * qsr(:,:) 
    134          ! 
    135          CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 )  
     136         ! 
     137         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
     138         ! 
     139         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )  
    136140         ! 
    137141         DO jk = 1, nksrp       
     
    161165         DO jj = 1, jpj 
    162166           DO ji = 1, jpi 
    163               IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN 
     167              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= zqsr100(ji,jj) )  THEN 
    164168                 neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer 
    165169                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint 
     
    226230      ENDIF 
    227231      ! 
    228       CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     232      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 ) 
     233      CALL wrk_dealloc( jpi, jpj,      zqsr100, zqsr_corr ) 
    229234      CALL wrk_dealloc( jpi, jpj, jpk, zpar,  ze0, ze1, ze2, ze3 ) 
    230235      ! 
     
    233238   END SUBROUTINE p4z_opt 
    234239 
    235    SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 )  
     240   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0, pqsr100 )  
    236241      !!---------------------------------------------------------------------- 
    237242      !!                  ***  routine p4z_opt_par  *** 
     
    242247      !!---------------------------------------------------------------------- 
    243248      !! * arguments 
    244       INTEGER, INTENT(in)                                       ::  kt            !   ocean time-step 
    245       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)              ::  pqsr          !   shortwave 
    246       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B) 
    247       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0   
     249      INTEGER, INTENT(in)                                        ::  kt            !   ocean time-step 
     250      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)               ::  pqsr          !   shortwave 
     251      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)            ::  pe1 , pe2 , pe3   !  PAR ( R-G-B) 
     252      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL  ::  pe0  
     253      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(out)  , OPTIONAL  ::  pqsr100   
    248254      !! * local variables 
    249255      INTEGER    ::   ji, jj, jk     ! dummy loop indices 
     
    255261      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:) 
    256262      ENDIF 
     263 
     264      !  Light at the euphotic depth  
     265      IF( PRESENT( pqsr100 ) )  pqsr100(:,:) = 0.01 * 3. * zqsr(:,:) 
    257266      ! 
    258267      IF( PRESENT( pe0 ) ) THEN     !  W-level 
     
    439448 
    440449   !!====================================================================== 
    441 END MODULE  p4zopt 
     450END MODULE p4zopt 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90

    r7959 r9987  
    202202                      zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj) 
    203203                      ! 
    204                       zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -znanotot ) ) 
    205                       zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trb(ji,jj,jk,jpdia) + rtrn ) 
    206  
    207                       zpislopen =  zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch)                & 
    208                         &          / ( trb(ji,jj,jk,jpphy) * 12.                  + rtrn )   & 
    209                         &          / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn ) 
    210  
    211                       zpislope2n = zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch)                & 
    212                         &          / ( trb(ji,jj,jk,jpdia) * 12.                  + rtrn )   & 
    213                         &          / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn ) 
     204                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -znanotot ) )           & 
     205                         &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn) 
     206                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trb(ji,jj,jk,jpdia) + rtrn )   & 
     207                         &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn) 
    214208 
    215209                      ! Computation of production function for Carbon 
    216210                      !  --------------------------------------------- 
     211                      zpislopen  =  zpislopead(ji,jj,jk)  / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn ) 
     212                      zpislope2n =  zpislopead2(ji,jj,jk) / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn ) 
    217213                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) ) 
    218214                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) ) 
     
    220216                      !  Computation of production function for Chlorophyll 
    221217                      !-------------------------------------------------- 
    222                       zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) ) 
    223                       zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) ) 
     218                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) ) 
     219                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) ) 
    224220                  ENDIF 
    225221               END DO 
     
    227223         END DO 
    228224      ENDIF 
    229  
    230  
     225       
    231226      !  Computation of a proxy of the N/C ratio 
    232227      !  --------------------------------------- 
     
    278273            zmxltst = MAX( 0.e0, hmld(ji,jj) - heup(ji,jj) ) 
    279274            zmxlday = zmxltst * zmxltst * r1_rday 
    280             zmixnano(ji,jj) = 1. - zmxlday / ( 2. + zmxlday ) 
    281             zmixdiat(ji,jj) = 1. - zmxlday / ( 4. + zmxlday ) 
     275            zmixnano(ji,jj) = 1. - zmxlday / ( 1. + zmxlday ) 
     276            zmixdiat(ji,jj) = 1. - zmxlday / ( 2. + zmxlday ) 
    282277         END DO 
    283278      END DO 
    284279  
    285       !  Mixed-layer effect on production                                                                                
     280      !  Mixed-layer effect on production  
     281      !  Sea-ice effect on production 
     282 
    286283      DO jk = 1, jpkm1 
    287284         DO jj = 1, jpj 
     
    291288                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj) 
    292289               ENDIF 
     290                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
     291                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
    293292            END DO 
    294293         END DO 
     
    330329      END DO 
    331330 
    332       IF( ln_newprod ) THEN 
    333 !CDIR NOVERRCHK 
    334          DO jk = 1, jpkm1 
    335 !CDIR NOVERRCHK 
    336             DO jj = 1, jpj 
    337 !CDIR NOVERRCHK 
    338                DO ji = 1, jpi 
    339                   IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
    340                      zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj) 
    341                      zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj) 
    342                   ENDIF 
    343                   IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    344                      !  production terms for nanophyto. ( chlorophyll ) 
    345                      znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 
    346                      zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 
    347                      zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
    348                      zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / & 
    349                                         & (  zpislopead(ji,jj,jk) * znanotot +rtrn) 
    350                      !  production terms for diatomees ( chlorophyll ) 
    351                      zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 
    352                      zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 
    353                      zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
    354                      zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / & 
    355                                         & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn ) 
    356                   ENDIF 
    357                END DO 
    358             END DO 
    359          END DO 
    360       ELSE 
    361 !CDIR NOVERRCHK 
    362          DO jk = 1, jpkm1 
    363 !CDIR NOVERRCHK 
    364             DO jj = 1, jpj 
    365 !CDIR NOVERRCHK 
    366                DO ji = 1, jpi 
    367                   IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    368                      !  production terms for nanophyto. ( chlorophyll ) 
    369                      znanotot = enano(ji,jj,jk) 
    370                      zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trb(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk) 
    371                      zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
    372                      zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 144. * zprod            & 
    373                      &                    / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * znanotot +rtrn ) 
    374                      !  production terms for diatomees ( chlorophyll ) 
    375                      zdiattot = ediat(ji,jj,jk) 
    376                      zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trb(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk) 
    377                      zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
    378                      zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 144. * zprod             & 
    379                      &                    / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn ) 
    380                   ENDIF 
    381                END DO 
    382             END DO 
    383          END DO 
    384       ENDIF 
     331!CDIR NOVERRCHK 
     332      DO jk = 1, jpkm1 
     333!CDIR NOVERRCHK 
     334         DO jj = 1, jpj 
     335!CDIR NOVERRCHK 
     336            DO ji = 1, jpi 
     337               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
     338                  zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj) 
     339                  zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj) 
     340               ENDIF 
     341               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
     342                  !  production terms for nanophyto. ( chlorophyll ) 
     343                  znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 
     344                  zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 
     345                  zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
     346                  zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / & 
     347                                     & (  zpislopead(ji,jj,jk) * znanotot +rtrn) 
     348                  !  production terms for diatomees ( chlorophyll ) 
     349                  zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 
     350                  zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 
     351                  zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
     352                  zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / & 
     353                                     & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn ) 
     354               ENDIF 
     355            END DO 
     356         END DO 
     357      END DO 
    385358 
    386359      !   Update the arrays TRA which contain the biological sources and sinks 
     
    629602 
    630603   !!====================================================================== 
    631 END MODULE  p4zprod 
     604END MODULE p4zprod 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90

    r7959 r9987  
    4444   REAL(wp), PUBLIC ::  xsiremlab  !: fast remineralisation rate of POC  
    4545   REAL(wp), PUBLIC ::  xsilab     !: fraction of labile biogenic silica  
    46    REAL(wp), PUBLIC ::  oxymin     !: halk saturation constant for anoxia  
    47  
    4846 
    4947   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitr     !: denitrification array 
     
    111109                  zdepprod(ji,jj,jk) = zdepmin**0.273 
    112110               ENDIF 
    113             END DO 
    114          END DO 
    115       END DO 
    116  
    117       DO jk = 1, jpkm1 
    118          DO jj = 1, jpj 
    119             DO ji = 1, jpi 
    120                ! denitrification factor computed from O2 levels 
    121                nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    & 
    122                   &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  ) 
    123                nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) ) 
    124111            END DO 
    125112         END DO 
     
    357344      !! 
    358345      !!---------------------------------------------------------------------- 
    359       NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab,   & 
    360       &                   oxymin 
     346      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab 
    361347      INTEGER :: ios                 ! Local integer output status for namelist read 
    362348 
     
    380366         WRITE(numout,*) '    fraction of labile biogenic silica        xsilab    =', xsilab 
    381367         WRITE(numout,*) '    NH4 nitrification rate                    nitrif    =', nitrif 
    382          WRITE(numout,*) '    halk saturation constant for anoxia       oxymin    =', oxymin 
    383368      ENDIF 
    384369      ! 
    385       nitrfac (:,:,:) = 0._wp 
    386370      denitr  (:,:,:) = 0._wp 
    387371      denitnh4(:,:,:) = 0._wp 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90

    r7960 r9987  
    159159      IF( ln_ndepo ) THEN 
    160160         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN 
    161             CALL fld_read( kt, 1, sf_ndepo ) 
    162             DO jj = 1, jpj 
    163                DO ji = 1, jpi 
    164                   nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / rno3 / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn ) 
    165                END DO 
    166             END DO 
     161             zcoef = rno3 * 14E6 * ryyss 
     162             CALL fld_read( kt, 1, sf_ndepo ) 
     163             nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / fse3t(:,:,1)  
     164         ENDIF 
     165         IF( lk_vvl ) THEN 
     166           zcoef = rno3 * 14E6 * ryyss 
     167           nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / fse3t(:,:,1)  
    167168         ENDIF 
    168169      ENDIF 
     
    455456            DO jj = 1, jpj 
    456457               DO ji = 1, jpi 
    457                   zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) ) 
     458                  zexpide   = MIN( 8.,( gdept_0(ji,jj,jk) / 500. )**(-1.5) ) 
    458459                  zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2 
    459460                  zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 ) 
     
    465466         ironsed(:,:,jpk) = 0._wp 
    466467         DO jk = 1, jpkm1 
    467             ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday ) 
     468            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday ) 
    468469         END DO 
    469470         DEALLOCATE( zcmask) 
     
    483484         CALL iom_close( numhydro ) 
    484485         ! 
    485          hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp 
     486         DO jk = 1, jpk 
     487            hydrofe(:,:,jk) = ( hydrofe(:,:,jk) * hratio ) / ( e1e2t(:,:) * e3t_0(:,:,jk) * ryyss + rtrn ) / 1000._wp 
     488         ENDDO 
    486489         ! 
    487490      ENDIF 
     
    519522 
    520523   !!====================================================================== 
    521 END MODULE  p4zsbc 
     524END MODULE p4zsbc 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90

    r7960 r9987  
    7272      CHARACTER (len=25) :: charout 
    7373      REAL(wp), POINTER, DIMENSION(:,:  ) :: zpdep, zsidep, zwork1, zwork2, zwork3 
     74      REAL(wp), POINTER, DIMENSION(:,:)   :: zsedcal, zsedsi, zsedc 
    7475      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdenit2d, zironice, zbureff 
    7576      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsbio3, zwsbio4, zwscal 
     
    8384      ! Allocate temporary workspace 
    8485      CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
     86      CALL wrk_alloc( jpi, jpj, zsedcal,  zsedsi, zsedc ) 
    8587      CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
    8688      CALL wrk_alloc( jpi, jpj, jpk, zsoufer ) 
     
    9193      zwork2  (:,:) = 0.e0 
    9294      zwork3  (:,:) = 0.e0 
     95      zsedsi   (:,:) = 0.e0 
     96      zsedcal  (:,:) = 0.e0 
     97      zsedc    (:,:) = 0.e0 
    9398 
    9499      ! Iron input/uptake due to sea ice : Crude parameterization based on Lancelot et al. 
     
    298303            tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
    299304            tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
     305            zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss / zdep 
     306            zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss / zdep 
    300307#endif 
    301308         END DO 
     
    336343            tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) ) 
    337344            tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt 
    338             sdenit(ji,jj) = rdenit * zpdenit * fse3t(ji,jj,ikt) 
     345            sdenit(ji,jj) = rdenit * zpdenit / zdep 
     346            zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc / zdep 
    339347#endif 
    340348         END DO 
     
    392400               CALL iom_put( "INTNFIX" , zwork1 )  
    393401            ENDIF 
     402            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * 1.e+3 ) 
     403            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * 1.e+3 ) 
     404            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * 1.e+3 ) 
     405            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * 1.e+3 * rno3 ) 
    394406         ENDIF 
    395407      ELSE 
     
    405417      ! 
    406418      CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
     419      CALL wrk_dealloc( jpi, jpj, zsedcal , zsedsi, zsedc ) 
    407420      CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
    408421      CALL wrk_dealloc( jpi, jpj, jpk, zsoufer ) 
     
    436449 
    437450   !!====================================================================== 
    438 END MODULE  p4zsed 
     451END MODULE p4zsed 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90

    r7959 r9987  
    913913 
    914914   !!====================================================================== 
    915 END MODULE  p4zsink 
     915END MODULE p4zsink 
  • branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90

    r7959 r9987  
    3838 
    3939   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget, po4budget 
    40    REAL(wp) :: xfact1, xfact2 
     40   REAL(wp) :: xfact1, xfact2, xfact3 
    4141   INTEGER ::  numco2, numnut, numnit  !: logical unit for co2 budget 
    4242 
     
    133133         ! 
    134134         CALL p4z_bio( kt, jnt )   ! Biology 
    135          CALL p4z_sed( kt, jnt )   ! Sedimentation 
    136135         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation 
     136         CALL p4z_sed( kt, jnt )   ! Surface and Bottom boundary conditions 
    137137         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes 
    138138         ! 
     
    474474      !!--------------------------------------------------------------------- 
    475475      ! 
    476       INTEGER , INTENT( in ) ::   kt      ! ocean time-step index       
    477       REAL(wp)               ::  zfact        
    478       REAL(wp) ::  zrdenittot, zsdenittot, znitrpottot 
     476      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index       
     477      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot 
    479478      CHARACTER(LEN=100)   ::   cltxt 
    480479      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol 
     
    492491            xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr 
    493492            xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr 
     493            xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s 
    494494            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron' 
    495495            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt) 
     
    574574      IF( iom_use( "Sdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN 
    575575         zsdenittot   = glob_sum ( sdenit(:,:) * e1e2t(:,:) ) 
    576          CALL iom_put( "Sdenit", sdenit(:,:) * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments 
     576         CALL iom_put( "Sdenit", sdenit(:,:) * xfact3 * tmask(:,:,1) )  ! Nitrate reduction in the sediments 
    577577      ENDIF 
    578578 
Note: See TracChangeset for help on using the changeset viewer.