Changeset 6945 for trunk


Ignore:
Timestamp:
2016-09-23T12:31:28+02:00 (5 years ago)
Author:
cetlod
Message:

trunk: bugfixes on PISCES carbon chemistry, see ticket #1774

Location:
trunk/NEMOGCM
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/CONFIG/SHARED/namelist_pisces_ref

    r5385 r6945  
    6666   qdfelim    =  7.E-6    ! Optimal quota of diatoms 
    6767   caco3r     =  0.3      ! mean rain ratio 
     68   oxymin    =  1.E-6     ! Half-saturation constant for anoxia 
    6869/ 
    6970!''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 
     
    162163   xsiremlab =  0.03      ! fast remineralization rate of Si 
    163164   xsilab    =  0.5       ! Fraction of labile biogenic silica 
    164    oxymin    =  1.E-6     ! Half-saturation constant for anoxia 
    165165/ 
    166166!''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90

    r6291 r6945  
    3131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sio3eq   ! chemistry of Si 
    3232   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fekeq    ! chemistry of Fe 
    33    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ::   chemc    ! Solubilities of O2 and CO2 
     33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemc    ! Solubilities of O2 and CO2 
    3434   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 
    5445   REAL(wp) ::   bor1   = 0.00023        ! borat constants 
    5546   REAL(wp) ::   bor2   = 1. / 10.82 
    56  
    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 
    7547 
    7648   REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate 
     
    144116      REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 
    145117      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2 
    146       REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1 
    147       REAL(wp) ::   zis  , zis2 , zsal15, zisqrt 
     118      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat 
     119      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1  , za2 
    148120      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw 
    149121      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1 
     
    151123      ! 
    152124      IF( nn_timing == 1 )  CALL timing_start('p4z_che') 
     125      ! 
     126      ! Computations of chemical constants require in situ temperature 
     127      ! Here a quite simple formulation is used to convert  
     128      ! potential temperature to in situ temperature. The errors is less than  
     129      ! 0.04°C relative to an exact computation 
     130      ! --------------------------------------------------------------------- 
     131      DO jk = 1, jpk 
     132         DO jj = 1, jpj 
     133            DO ji = 1, jpi 
     134               zpres = gdept_n(ji,jj,jk) / 1000. 
     135               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (tsn(ji,jj,jk,jp_sal) - 35.0) ) 
     136               za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 ) 
     137               tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2 
     138            END DO 
     139         END DO 
     140      END DO 
    153141      ! 
    154142      ! CHEMICAL CONSTANTS - SURFACE LAYER 
     
    157145         DO ji = 1, jpi 
    158146            !                             ! SET ABSOLUTE TEMPERATURE 
    159             ztkel = tsn(ji,jj,1,jp_tem) + 273.15 
     147            ztkel = tempis(ji,jj,1) + 273.15 
    160148            zt    = ztkel * 0.01 
    161149            zt2   = zt * zt 
     
    165153            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 
    166154            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
    167             zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 ) 
     155            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    & 
     156            &       + 0.0047036e-4*ztkel**2) 
    168157            !                             ! SET SOLUBILITIES OF O2 AND CO2  
    169             chemc(ji,jj) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm) 
     158            chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm) 
     159            chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 
     160            chemc(ji,jj,3) = 57.7 - 0.118*ztkel 
    170161            ! 
    171162         END DO 
     
    177168         DO jj = 1, jpj 
    178169            DO ji = 1, jpi 
    179               ztkel = tsn(ji,jj,jk,jp_tem) + 273.15 
     170              ztkel = tempis(ji,jj,jk) + 273.15 
    180171              zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35. 
    181172              zsal2 = zsal * zsal 
    182               ztgg  = LOG( ( 298.15 - tsn(ji,jj,jk,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
     173              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
    183174              ztgg2 = ztgg  * ztgg 
    184175              ztgg3 = ztgg2 * ztgg 
     
    200191            DO ji = 1, jpi 
    201192 
    202                ! SET PRESSION 
    203                zpres   = 1.025e-1 * gdept_n(ji,jj,jk) 
     193               ! SET PRESSION ACCORDING TO SAUNDER (1980) 
     194               zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) ) 
     195               zc1 = 5.92E-3 + zplat**2 * 5.25E-3 
     196               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept_n(ji,jj,jk)))) / 4.42E-6 
     197               zpres = zpres / 10.0 
    204198 
    205199               ! SET ABSOLUTE TEMPERATURE 
    206                ztkel   = tsn(ji,jj,jk,jp_tem) + 273.15 
     200               ztkel   = tempis(ji,jj,jk) + 273.15 
    207201               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35. 
    208202               zsqrt  = SQRT( zsal ) 
     
    213207               zis2   = zis * zis 
    214208               zisqrt = SQRT( zis ) 
    215                ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20. 
     209               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20. 
    216210 
    217211               ! CHLORINITY (WOOSTER ET AL., 1969) 
     
    246240 
    247241 
    248                zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal 
    249                zck2    = c20 * ztr + c21 + c22 * zsal   + c23 * zsal**2 
     242               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO  
     243               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 
     244               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  & 
     245                  - 0.011555*zsal + 0.0001152*zsal*zsal) 
     246               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      & 
     247                  - 0.01781*zsal + 0.0001122*zsal*zsal) 
    250248 
    251249               ! PKW (H2O) (DICKSON AND RILEY, 1979) 
     
    256254               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
    257255               !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 
    258                zaksp0  = akcc1 + akcc2 * ztkel + akcc3 * ztr + akcc4 * LOG10( ztkel )   & 
    259                   &   + ( akcc5 + akcc6 * ztkel + akcc7 * ztr ) * zsqrt + akcc8 * zsal + akcc9 * zsal15 
     256               zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   & 
     257                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  & 
     258                  &      - 0.07711*zsal + 0.0041249*zsal15 
    260259 
    261260               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 
     
    327326      !!                     ***  ROUTINE p4z_che_alloc  *** 
    328327      !!---------------------------------------------------------------------- 
    329       ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj), chemo2(jpi,jpj,jpk),   & 
    330       &         STAT=p4z_che_alloc ) 
     328      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk),   & 
     329      &         tempis(jpi,jpj,jpk), STAT=p4z_che_alloc ) 
    331330      ! 
    332331      IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90

    r6291 r6945  
    8484      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan 
    8585      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact 
     86      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff 
    8687      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 
    8788      REAL(wp) ::   zyr_dec, zdco2dt 
    8889      CHARACTER (len=25) :: charout 
    89       REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d  
     90      REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx, zw2d, zpco2atm  
    9091      !!--------------------------------------------------------------------- 
    9192      ! 
    9293      IF( nn_timing == 1 )  CALL timing_start('p4z_flx') 
    9394      ! 
    94       CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx ) 
     95      CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 
    9596      ! 
    9697 
     
    177178      DO jj = 1, jpj 
    178179         DO ji = 1, jpi 
     180            ztkel     = tsn(ji,jj,1,jp_tem) + 273.15 
     181            zsal      = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35. 
     182            zvapsw    = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal) 
     183            zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw ) 
     184            zxc2      = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2 
     185            zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) )   & 
     186            &           / ( 82.05736 * ztkel )) 
     187            zfco2 = zpco2atm(ji,jj) * zfugcoeff 
     188 
    179189            ! Compute CO2 flux for the sea and air 
    180             zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj) * zkgco2(ji,jj)   ! (mol/L) * (m/s) 
    181             zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ? 
     190            zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj)  ! (mol/L) * (m/s) 
     191            zflu = zh2co3(ji,jj) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ? 
    182192            oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 
    183193            ! compute the trend 
    184             tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / e3t_n(ji,jj,1) 
     194            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / e3t_n(ji,jj,1) * tmask(ji,jj,1) 
    185195 
    186196            ! Compute O2 flux  
    187             zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
    188             zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj) 
    189             zoflx(ji,jj) = zfld16 - zflu16 
     197            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
     198            zflu16 = trb(ji,jj,1,jpoxy) * zkgo2(ji,jj) 
     199            zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1) 
    190200            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / e3t_n(ji,jj,1) 
    191201         END DO 
     
    218228         ENDIF 
    219229         IF( iom_use( "Dpco2" ) ) THEN 
    220            zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1) 
     230           zw2d(:,:) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1) 
    221231           CALL iom_put( "Dpco2" ,  zw2d ) 
    222232         ENDIF 
     
    234244            trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1)  
    235245            trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1)  
    236             trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1)  
    237          ENDIF 
    238       ENDIF 
    239       ! 
    240       CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx ) 
     246            trc2d(:,:,jp_pcs0_2d + 3) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1)  
     247         ENDIF 
     248      ENDIF 
     249      ! 
     250      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 
    241251      ! 
    242252      IF( nn_timing == 1 )  CALL timing_stop('p4z_flx') 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90

    r5836 r6945  
    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 
     
    186187      END DO 
    187188      ! 
     189      DO jk = 1, jpkm1 
     190         DO jj = 1, jpj 
     191            DO ji = 1, jpi 
     192               ! denitrification factor computed from O2 levels 
     193               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    & 
     194                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  ) 
     195               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) ) 
     196            END DO 
     197         END DO 
     198      END DO 
    188199      ! 
    189200      IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics 
     
    215226      NAMELIST/nampislim/ concnno3, concdno3, concnnh4, concdnh4, concnfer, concdfer, concbfe,   & 
    216227         &                concbno3, concbnh4, xsizedia, xsizephy, xsizern, xsizerd,          &  
    217          &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r 
     228         &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r, oxymin 
    218229      INTEGER :: ios                 ! Local integer output status for namelist read 
    219230 
     
    248259         WRITE(numout,*) '    Minimum size criteria for nanophyto      xsizephy  = ', xsizephy 
    249260         WRITE(numout,*) '    Fe half saturation for bacteria          concbfe   = ', concbfe 
     261         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin   =' , oxymin 
    250262         WRITE(numout,*) '    optimal Fe quota for nano.               qnfelim   = ', qnfelim 
    251263         WRITE(numout,*) '    Optimal Fe quota for diatoms             qdfelim   = ', qdfelim 
    252264      ENDIF 
    253  
     265      ! 
     266      nitrfac (:,:,:) = 0._wp 
     267      ! 
    254268   END SUBROUTINE p4z_lim_init 
    255269 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

    r6291 r6945  
    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. 
     
    117117               zcalcon  = calcon * ( tsn(ji,jj,jk,jp_sal) / 35._wp ) 
    118118               zfact    = rhop(ji,jj,jk) / 1000._wp 
    119                zomegaca = ( zcalcon * zco3(ji,jj,jk) * zfact ) / aksp(ji,jj,jk)  
     119               zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn ) 
     120               zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn ) 
    120121 
    121122               ! SET DEGREE OF UNDER-/SUPERSATURATION 
     
    146147      IF( lk_iomput .AND. knt == nrdttrc ) THEN 
    147148         IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )          * tmask(:,:,:) ) 
    148          IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:) * 1.e+3               * tmask(:,:,:) ) 
    149          IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", aksp(:,:,:) * 1.e+3 / calcon      * tmask(:,:,:) ) 
    150          IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) ) 
     149         IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:)    * 1.e+3            * tmask(:,:,:) ) 
     150         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", zco3sat(:,:,:) * 1.e+3            * tmask(:,:,:) ) 
     151         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) ) 
    151152      ELSE 
    152153         IF( ln_diatrc ) THEN 
    153154            trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:) 
    154155            trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:) 
    155             trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:) 
     156            trc3d(:,:,:,jp_pcs0_3d + 2) = zco3sat(:,:,:)           * tmask(:,:,:) 
    156157         ENDIF 
    157158      ENDIF 
     
    163164      ENDIF 
    164165      ! 
    165       CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss ) 
     166      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss ) 
    166167      ! 
    167168      IF( nn_timing == 1 )  CALL timing_stop('p4z_lys') 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90

    r6628 r6945  
    110110      !                                        !  -------------------------------------- 
    111111      IF( l_trcdm2dc ) THEN                     !  diurnal cycle 
    112          !                                       ! 1% of qsr to compute euphotic layer 
    113          zqsr100(:,:)   = 0.01 * qsr_mean(:,:)     !  daily mean qsr 
    114112         ! 
    115113         zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
    116114         ! 
    117          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )  
     115         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )  
    118116         ! 
    119117         DO jk = 1, nksrp       
     
    132130         ! 
    133131      ELSE 
    134          !                                       ! 1% of qsr to compute euphotic layer 
    135          zqsr100(:,:)   = 0.01 * qsr(:,:)     !  daily mean qsr 
    136132         ! 
    137133         zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn ) 
    138134         ! 
    139          CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )  
     135         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 )  
    140136         ! 
    141137         DO jk = 1, nksrp       
     
    165161         DO jj = 1, jpj 
    166162           DO ji = 1, jpi 
    167               IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN 
     163              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= zqsr100(ji,jj) )  THEN 
    168164                 neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer 
    169165                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint 
     
    234230   END SUBROUTINE p4z_opt 
    235231 
    236    SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 )  
     232   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0, pqsr100 )  
    237233      !!---------------------------------------------------------------------- 
    238234      !!                  ***  routine p4z_opt_par  *** 
     
    247243      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B) 
    248244      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0   
     245      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out)  , OPTIONAL  ::  pqsr100   
    249246      !! * local variables 
    250247      INTEGER    ::   ji, jj, jk     ! dummy loop indices 
     
    256253      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:) 
    257254      ENDIF 
    258       ! 
     255       
     256      !  Light at the euphotic depth  
     257      IF( PRESENT( pqsr100 ) )  pqsr100(:,:) = 0.01 * 3. * zqsr(:,:) 
     258 
    259259      IF( PRESENT( pe0 ) ) THEN     !  W-level 
    260260         ! 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90

    r6628 r6945  
    279279                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj) 
    280280               ENDIF 
     281               zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
     282               zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
    281283            END DO 
    282284         END DO 
     
    315317      END DO 
    316318 
    317       IF( ln_newprod ) THEN 
    318          DO jk = 1, jpkm1 
    319             DO jj = 1, jpj 
    320                DO ji = 1, jpi 
    321                   IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
    322                      zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj) 
    323                      zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj) 
    324                   ENDIF 
    325                   IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    326                      !  production terms for nanophyto. ( chlorophyll ) 
    327                      znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 
    328                      zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 
    329                      zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
    330                      zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / & 
    331                                         & (  zpislopead(ji,jj,jk) * znanotot +rtrn) 
    332                      !  production terms for diatomees ( chlorophyll ) 
    333                      zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 
    334                      zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 
    335                      zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
    336                      zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / & 
    337                                         & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn ) 
    338                   ENDIF 
    339                END DO 
    340             END DO 
    341          END DO 
    342       ELSE 
    343          DO jk = 1, jpkm1 
    344             DO jj = 1, jpj 
    345                DO ji = 1, jpi 
    346                   IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
    347                      !  production terms for nanophyto. ( chlorophyll ) 
    348                      znanotot = enano(ji,jj,jk) 
    349                      zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trb(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk) 
    350                      zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
    351                      zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 144. * zprod            & 
    352                      &                    / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * znanotot +rtrn ) 
    353                      !  production terms for diatomees ( chlorophyll ) 
    354                      zdiattot = ediat(ji,jj,jk) 
    355                      zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trb(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk) 
    356                      zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
    357                      zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 144. * zprod             & 
    358                      &                    / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn ) 
    359                   ENDIF 
    360                END DO 
    361             END DO 
    362          END DO 
    363       ENDIF 
     319      DO jk = 1, jpkm1 
     320         DO jj = 1, jpj 
     321            DO ji = 1, jpi 
     322               IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN 
     323                  zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj) 
     324                  zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj) 
     325               ENDIF 
     326               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN 
     327                  !  production terms for nanophyto. ( chlorophyll ) 
     328                  znanotot = enano(ji,jj,jk) * zstrn(ji,jj) 
     329                  zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk) 
     330                  zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk) 
     331                  zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / & 
     332                                     & (  zpislopead(ji,jj,jk) * znanotot +rtrn) 
     333                  !  production terms for diatomees ( chlorophyll ) 
     334                  zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj) 
     335                  zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk) 
     336                  zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk) 
     337                  zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / & 
     338                                     & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn ) 
     339               ENDIF 
     340            END DO 
     341         END DO 
     342      END DO 
    364343 
    365344      !   Update the arrays TRA which contain the biological sources and sinks 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90

    r6140 r6945  
    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  
    4746 
    4847 
     
    109108                  zdepprod(ji,jj,jk) = zdepmin**0.273 
    110109               ENDIF 
    111             END DO 
    112          END DO 
    113       END DO 
    114  
    115       DO jk = 1, jpkm1 
    116          DO jj = 1, jpj 
    117             DO ji = 1, jpi 
    118                ! denitrification factor computed from O2 levels 
    119                nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    & 
    120                   &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  ) 
    121                nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) ) 
    122110            END DO 
    123111         END DO 
     
    355343      !! 
    356344      !!---------------------------------------------------------------------- 
    357       NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab,   & 
    358       &                   oxymin 
     345      NAMELIST/nampisrem/ xremik, xremip, nitrif, xsirem, xsiremlab, xsilab 
    359346      INTEGER :: ios                 ! Local integer output status for namelist read 
    360347 
     
    378365         WRITE(numout,*) '    fraction of labile biogenic silica        xsilab    =', xsilab 
    379366         WRITE(numout,*) '    NH4 nitrification rate                    nitrif    =', nitrif 
    380          WRITE(numout,*) '    halk saturation constant for anoxia       oxymin    =', oxymin 
    381367      ENDIF 
    382368      ! 
    383       nitrfac (:,:,:) = 0._wp 
    384369      denitr  (:,:,:) = 0._wp 
    385370      denitnh4(:,:,:) = 0._wp 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90

    r6140 r6945  
    155155      IF( ln_ndepo ) THEN 
    156156         IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_ndep > 1 ) ) THEN 
    157             CALL fld_read( kt, 1, sf_ndepo ) 
    158             DO jj = 1, jpj 
    159                DO ji = 1, jpi 
    160                   nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / rno3 / ( 14E6 * ryyss * e3t_n(ji,jj,1) + rtrn ) 
    161                END DO 
    162             END DO 
     157             zcoef = rno3 * 14E6 * ryyss 
     158             CALL fld_read( kt, 1, sf_ndepo ) 
     159             nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / e3t_n(:,:,1)  
     160         ENDIF 
     161         IF( lk_vvl ) THEN 
     162           zcoef = rno3 * 14E6 * ryyss 
     163           nitdep(:,:) = sf_ndepo(1)%fnow(:,:,1) / zcoef / e3t_n(:,:,1)  
    163164         ENDIF 
    164165      ENDIF 
     
    461462         ironsed(:,:,jpk) = 0._wp 
    462463         DO jk = 1, jpkm1 
    463             ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_n(:,:,jk) * rday ) 
     464            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_0(:,:,jk) * rday ) 
    464465         END DO 
    465466         DEALLOCATE( zcmask) 
     
    479480         CALL iom_close( numhydro ) 
    480481         ! 
    481          hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp 
     482         DO jk = 1, jpk 
     483            hydrofe(:,:,jk) = ( hydrofe(:,:,jk) * hratio ) / ( e1e2t(:,:) * e3t_0(:,:,jk) * ryyss + rtrn ) / 1000._wp 
     484         ENDDO 
    482485         ! 
    483486      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.