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 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90 – NEMO

Ignore:
Timestamp:
2015-12-01T16:35:30+01:00 (8 years ago)
Author:
timgraham
Message:

Upgraded branch to r5518 of trunk (v3.6 stable revision)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90

    r4624 r5965  
    4141#endif 
    4242 
    43    INTEGER  :: iksed  = 10 
     43   INTEGER  :: ik100 
    4444 
    4545#if  defined key_kriest 
     
    7979   !!---------------------------------------------------------------------- 
    8080 
    81    SUBROUTINE p4z_sink ( kt, jnt ) 
     81   SUBROUTINE p4z_sink ( kt, knt ) 
    8282      !!--------------------------------------------------------------------- 
    8383      !!                     ***  ROUTINE p4z_sink  *** 
     
    8888      !! ** Method  : - ??? 
    8989      !!--------------------------------------------------------------------- 
    90       INTEGER, INTENT(in) :: kt, jnt 
     90      INTEGER, INTENT(in) :: kt, knt 
    9191      INTEGER  ::   ji, jj, jk, jit 
    9292      INTEGER  ::   iiter1, iiter2 
     
    9494      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3 
    9595      REAL(wp) ::   zfact, zwsmax, zmax, zstep 
    96       REAL(wp) ::   zrfact2 
    97       INTEGER  ::   ik1 
    9896      CHARACTER (len=25) :: charout 
     97      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d 
     98      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d 
    9999      !!--------------------------------------------------------------------- 
    100100      ! 
     
    199199               zfact = zstep * xdiss(ji,jj,jk) 
    200200               !  Part I : Coagulation dependent on turbulence 
    201                zagg1 = 25.9  * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) 
    202                zagg2 = 4452. * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) 
     201               zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc) 
     202               zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc) 
    203203 
    204204               ! Part II : Differential settling 
    205205 
    206206               !  Aggregation of small into large particles 
    207                zagg3 =  47.1 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) 
    208                zagg4 =  3.3  * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) 
     207               zagg3 =  47.1 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc) 
     208               zagg4 =  3.3  * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc) 
    209209 
    210210               zagg   = zagg1 + zagg2 + zagg3 + zagg4 
    211                zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn ) 
     211               zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn ) 
    212212 
    213213               ! Aggregation of DOC to POC :  
     
    215215               ! 2nd term is shear aggregation of DOC-POC 
    216216               ! 3rd term is differential settling of DOC-POC 
    217                zaggdoc  = ( ( 0.369 * 0.3 * trn(ji,jj,jk,jpdoc) + 102.4 * trn(ji,jj,jk,jppoc) ) * zfact       & 
    218                &            + 2.4 * zstep * trn(ji,jj,jk,jppoc) ) * 0.3 * trn(ji,jj,jk,jpdoc) 
     217               zaggdoc  = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       & 
     218               &            + 2.4 * zstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc) 
    219219               ! transfer of DOC to GOC :  
    220220               ! 1st term is shear aggregation 
    221221               ! 2nd term is differential settling  
    222                zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trn(ji,jj,jk,jpgoc) * 0.3 * trn(ji,jj,jk,jpdoc) 
     222               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc) 
    223223               ! tranfer of DOC to POC due to brownian motion 
    224                zaggdoc3 =  ( 5095. * trn(ji,jj,jk,jppoc) + 114. * 0.3 * trn(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trn(ji,jj,jk,jpdoc) 
     224               zaggdoc3 =  ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc) 
    225225 
    226226               !  Update the trends 
     
    235235      END DO 
    236236 
    237      ! Total primary production per year 
    238      t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,iksed+1) + sinking2(:,:,iksed+1) ) * e1e2t(:,:) * tmask(:,:,1) ) 
     237 
     238     ! Total carbon export per year 
     239     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  & 
     240        &   t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) ) 
    239241     ! 
    240      IF( ln_diatrc ) THEN 
    241          zrfact2 = 1.e3 * rfact2r 
    242          ik1  = iksed + 1 
    243          IF( lk_iomput ) THEN 
    244            IF( jnt == nrdttrc ) THEN 
    245               CALL iom_put( "EPC100"  , ( sinking(:,:,ik1) + sinking2(:,:,ik1) ) * zrfact2 * tmask(:,:,1) ) ! Export of carbon at 100m 
    246               CALL iom_put( "EPFE100" , ( sinkfer(:,:,ik1) + sinkfer2(:,:,ik1) ) * zrfact2 * tmask(:,:,1) ) ! Export of iron at 100m 
    247               CALL iom_put( "EPCAL100",   sinkcal(:,:,ik1)                       * zrfact2 * tmask(:,:,1) ) ! Export of calcite  at 100m 
    248               CALL iom_put( "EPSI100" ,   sinksil(:,:,ik1)                       * zrfact2 * tmask(:,:,1) ) ! Export of biogenic silica at 100m 
    249            ENDIF 
    250          ELSE 
    251            trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik1) * zrfact2 * tmask(:,:,1) 
    252            trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik1) * zrfact2 * tmask(:,:,1) 
    253            trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik1) * zrfact2 * tmask(:,:,1) 
    254            trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik1) * zrfact2 * tmask(:,:,1) 
    255            trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik1) * zrfact2 * tmask(:,:,1) 
    256            trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik1) * zrfact2 * tmask(:,:,1) 
     242     IF( lk_iomput ) THEN 
     243       IF( knt == nrdttrc ) THEN 
     244          CALL wrk_alloc( jpi, jpj,      zw2d ) 
     245          CALL wrk_alloc( jpi, jpj, jpk, zw3d ) 
     246          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s 
     247          ! 
     248          IF( iom_use( "EPC100" ) )  THEN 
     249              zw2d(:,:) = ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of carbon at 100m 
     250              CALL iom_put( "EPC100"  , zw2d ) 
     251          ENDIF 
     252          IF( iom_use( "EPFE100" ) )  THEN 
     253              zw2d(:,:) = ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of iron at 100m 
     254              CALL iom_put( "EPFE100"  , zw2d ) 
     255          ENDIF 
     256          IF( iom_use( "EPCAL100" ) )  THEN 
     257              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m 
     258              CALL iom_put( "EPCAL100"  , zw2d ) 
     259          ENDIF 
     260          IF( iom_use( "EPSI100" ) )  THEN 
     261              zw2d(:,:) =  sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m 
     262              CALL iom_put( "EPSI100"  , zw2d ) 
     263          ENDIF 
     264          IF( iom_use( "EXPC" ) )  THEN 
     265              zw3d(:,:,:) = ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of carbon in the water column 
     266              CALL iom_put( "EXPC"  , zw3d ) 
     267          ENDIF 
     268          IF( iom_use( "EXPFE" ) )  THEN 
     269              zw3d(:,:,:) = ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of iron  
     270              CALL iom_put( "EXPFE"  , zw3d ) 
     271          ENDIF 
     272          IF( iom_use( "EXPCAL" ) )  THEN 
     273              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite  
     274              CALL iom_put( "EXPCAL"  , zw3d ) 
     275          ENDIF 
     276          IF( iom_use( "EXPSI" ) )  THEN 
     277              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica 
     278              CALL iom_put( "EXPSI"  , zw3d ) 
     279          ENDIF 
     280          IF( iom_use( "tcexp" ) )  CALL iom_put( "tcexp" , t_oce_co2_exp * zfact )   ! molC/s 
     281          !  
     282          CALL wrk_dealloc( jpi, jpj,      zw2d ) 
     283          CALL wrk_dealloc( jpi, jpj, jpk, zw3d ) 
     284        ENDIF 
     285      ELSE 
     286         IF( ln_diatrc ) THEN 
     287            zfact = 1.e3 * rfact2r 
     288            trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1) 
     289            trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1) 
     290            trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1) 
     291            trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1) 
     292            trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1) 
     293            trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1) 
    257294         ENDIF 
    258295      ENDIF 
     
    272309      !!                  ***  ROUTINE p4z_sink_init  *** 
    273310      !!---------------------------------------------------------------------- 
    274  
     311      INTEGER :: jk 
     312 
     313      ik100 = 10        !  last level where depth less than 100 m 
     314      DO jk = jpkm1, 1, -1 
     315         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1 
     316      END DO 
     317      IF (lwp) WRITE(numout,*) 
     318      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1 
     319      IF (lwp) WRITE(numout,*) 
     320      ! 
    275321      t_oce_co2_exp = 0._wp 
    276322      ! 
     
    282328   !!---------------------------------------------------------------------- 
    283329 
    284    SUBROUTINE p4z_sink ( kt, jnt ) 
     330   SUBROUTINE p4z_sink ( kt, knt ) 
    285331      !!--------------------------------------------------------------------- 
    286332      !!                ***  ROUTINE p4z_sink  *** 
     
    292338      !!--------------------------------------------------------------------- 
    293339      ! 
    294       INTEGER, INTENT(in) :: kt, jnt 
     340      INTEGER, INTENT(in) :: kt, knt 
    295341      ! 
    296342      INTEGER  :: ji, jj, jk, jit, niter1, niter2 
     
    300346      REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5 
    301347      REAL(wp) :: zval1, zval2, zval3, zval4 
    302       REAL(wp) :: zrfact2 
     348      REAL(wp) :: zfact 
    303349      INTEGER  :: ik1 
    304350      CHARACTER (len=25) :: charout 
    305351      REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d  
     352      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d 
     353      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d 
    306354      !!--------------------------------------------------------------------- 
    307355      ! 
     
    325373            DO ji = 1, jpi 
    326374               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN 
    327                   znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp 
     375                  znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp 
    328376                  ! -------------- To avoid sinking speed over 50 m/day ------- 
    329377                  znum  = MIN( xnumm(jk), znum ) 
     
    387435               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN 
    388436 
    389                   znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp 
     437                  znum = trb(ji,jj,jk,jppoc)/(trb(ji,jj,jk,jpnum)+rtrn) / xkr_massp 
    390438                  !-------------- To avoid sinking speed over 50 m/day ------- 
    391439                  znum  = min(xnumm(jk),znum) 
     
    405453                  !    ---------------------------------------------- 
    406454 
    407                   zagg1 =  0.163 * trn(ji,jj,jk,jpnum)**2               & 
     455                  zagg1 =  0.163 * trb(ji,jj,jk,jpnum)**2               & 
    408456                     &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    & 
    409457                     &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    & 
    410458                     &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  & 
    411459                     &            * (zeps-1.)**2/(zdiv2*zdiv3))  
    412                   zagg2 =  2*0.163*trn(ji,jj,jk,jpnum)**2*zfm*                       & 
     460                  zagg2 =  2*0.163*trb(ji,jj,jk,jpnum)**2*zfm*                       & 
    413461                     &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2          & 
    414462                     &                    *xkr_mass_min*(zeps-1.)/zdiv2                 & 
     
    418466                     &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))     
    419467 
    420                   zagg3 =  0.163*trn(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3   
     468                  zagg3 =  0.163*trb(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3   
    421469                   
    422470                 !    Aggregation of small into large particles 
     
    424472                 !    ---------------------------------------------- 
    425473 
    426                   zagg4 =  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2*                       & 
     474                  zagg4 =  2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2*                       & 
    427475                     &                 xkr_wsbio_min*(zeps-1.)**2                         & 
    428476                     &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      & 
     
    431479                     &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )    
    432480 
    433                   zagg5 =   2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2                         & 
     481                  zagg5 =   2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2                         & 
    434482                     &                 *(zeps-1.)*zfm*xkr_wsbio_min                        & 
    435483                     &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         & 
     
    441489                  !     ------------------------------------ 
    442490 
    443                   zfract = 2.*3.141*0.125*trn(ji,jj,jk,jpmes)*12./0.12/0.06**3*trn(ji,jj,jk,jpnum)  & 
     491                  zfract = 2.*3.141*0.125*trb(ji,jj,jk,jpmes)*12./0.12/0.06**3*trb(ji,jj,jk,jpnum)  & 
    444492                    &      * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2  & 
    445493                    &      * 10000.*xstep 
     
    448496                  !     -------------------------------------- 
    449497 
    450                   zaggdoc = 0.83 * trn(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)   & 
    451                      &        + 0.005 * 231. * trn(ji,jj,jk,jpdoc) * xstep * trn(ji,jj,jk,jpdoc) 
    452                   zaggdoc1 = 271. * trn(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)  & 
    453                      &  + 0.02 * 16706. * trn(ji,jj,jk,jppoc) * xstep * trn(ji,jj,jk,jpdoc) 
     498                  zaggdoc = 0.83 * trb(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)   & 
     499                     &        + 0.005 * 231. * trb(ji,jj,jk,jpdoc) * xstep * trb(ji,jj,jk,jpdoc) 
     500                  zaggdoc1 = 271. * trb(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)  & 
     501                     &  + 0.02 * 16706. * trb(ji,jj,jk,jppoc) * xstep * trb(ji,jj,jk,jpdoc) 
    454502 
    455503# if defined key_degrad 
     
    466514                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi ) 
    467515                  ! 
    468                   znumdoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn ) 
     516                  znumdoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn ) 
    469517                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1 
    470518                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg 
     
    477525 
    478526     ! Total primary production per year 
    479      t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,:) ) * cvol(:,:,:) ) 
     527     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) ) 
    480528     ! 
    481       IF( ln_diatrc ) THEN 
    482          ! 
    483          ik1 = iksed + 1 
    484          zrfact2 = 1.e3 * rfact2r 
    485          IF( jnt == nrdttrc ) THEN 
    486            CALL iom_put( "POCFlx"  , sinking (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! POC export 
    487            CALL iom_put( "NumFlx"  , sinking2 (:,:,:)     * zrfact2 * tmask(:,:,:) )  ! Num export 
    488            CALL iom_put( "SiFlx"   , sinksil (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! Silica export 
    489            CALL iom_put( "CaCO3Flx", sinkcal (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! Calcite export 
    490            CALL iom_put( "xnum"    , znum3d  (:,:,:)                * tmask(:,:,:) )  ! Number of particles in aggregats 
    491            CALL iom_put( "W1"      , wsbio3  (:,:,:)                * tmask(:,:,:) )  ! sinking speed of POC 
    492            CALL iom_put( "W2"      , wsbio4  (:,:,:)                * tmask(:,:,:) )  ! sinking speed of aggregats 
     529     IF( lk_iomput ) THEN 
     530        IF( knt == nrdttrc ) THEN 
     531          CALL wrk_alloc( jpi, jpj,      zw2d ) 
     532          CALL wrk_alloc( jpi, jpj, jpk, zw3d ) 
     533          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s 
     534          ! 
     535          IF( iom_use( "EPC100" ) )  THEN 
     536              zw2d(:,:) = sinking(:,:,ik100) * zfact * tmask(:,:,1) ! Export of carbon at 100m 
     537              CALL iom_put( "EPC100"  , zw2d ) 
     538          ENDIF 
     539          IF( iom_use( "EPN100" ) )  THEN 
     540              zw2d(:,:) = sinking2(:,:,ik100) * zfact * tmask(:,:,1) ! Export of number of aggregates ? 
     541              CALL iom_put( "EPN100"  , zw2d ) 
     542          ENDIF 
     543          IF( iom_use( "EPCAL100" ) )  THEN 
     544              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m 
     545              CALL iom_put( "EPCAL100"  , zw2d ) 
     546          ENDIF 
     547          IF( iom_use( "EPSI100" ) )  THEN 
     548              zw2d(:,:) = sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m 
     549              CALL iom_put( "EPSI100"  , zw2d ) 
     550          ENDIF 
     551          IF( iom_use( "EXPC" ) )  THEN 
     552              zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column 
     553              CALL iom_put( "EXPC"  , zw3d ) 
     554          ENDIF 
     555          IF( iom_use( "EXPN" ) )  THEN 
     556              zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column 
     557              CALL iom_put( "EXPN"  , zw3d ) 
     558          ENDIF 
     559          IF( iom_use( "EXPCAL" ) )  THEN 
     560              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite  
     561              CALL iom_put( "EXPCAL"  , zw3d ) 
     562          ENDIF 
     563          IF( iom_use( "EXPSI" ) )  THEN 
     564              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica 
     565              CALL iom_put( "EXPSI"  , zw3d ) 
     566          ENDIF 
     567          IF( iom_use( "XNUM" ) )  THEN 
     568              zw3d(:,:,:) =  znum3d(:,:,:) * tmask(:,:,:) !  Number of particles on aggregats 
     569              CALL iom_put( "XNUM"  , zw3d ) 
     570          ENDIF 
     571          IF( iom_use( "WSC" ) )  THEN 
     572              zw3d(:,:,:) = wsbio3(:,:,:) * tmask(:,:,:) ! Sinking speed of carbon particles 
     573              CALL iom_put( "WSC"  , zw3d ) 
     574          ENDIF 
     575          IF( iom_use( "WSN" ) )  THEN 
     576              zw3d(:,:,:) = wsbio4(:,:,:) * tmask(:,:,:) ! Sinking speed of particles number 
     577              CALL iom_put( "WSN"  , zw3d ) 
     578          ENDIF 
     579          ! 
     580          CALL wrk_dealloc( jpi, jpj,      zw2d ) 
     581          CALL wrk_dealloc( jpi, jpj, jpk, zw3d ) 
     582      ELSE 
     583         IF( ln_diatrc ) THEN 
     584            zfact = 1.e3 * rfact2r 
     585            trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,ik100)  * zfact * tmask(:,:,1) 
     586            trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,ik100)  * zfact * tmask(:,:,1) 
     587            trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,ik100)  * zfact * tmask(:,:,1) 
     588            trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,ik100)  * zfact * tmask(:,:,1) 
     589            trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,ik100)  * zfact * tmask(:,:,1) 
     590            trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:)      * zfact * tmask(:,:,:) 
     591            trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:)      * zfact * tmask(:,:,:) 
     592            trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:)      * zfact * tmask(:,:,:) 
     593            trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:)      * zfact * tmask(:,:,:) 
     594            trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d  (:,:,:)              * tmask(:,:,:) 
     595            trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3  (:,:,:)              * tmask(:,:,:) 
     596            trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4  (:,:,:)              * tmask(:,:,:) 
    493597         ENDIF 
    494 # if ! defined key_iomput 
    495          trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,ik1)    * zrfact2 * tmask(:,:,1) 
    496          trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,ik1)    * zrfact2 * tmask(:,:,1) 
    497          trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,ik1)    * zrfact2 * tmask(:,:,1) 
    498          trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,ik1)    * zrfact2 * tmask(:,:,1) 
    499          trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,ik1)    * zrfact2 * tmask(:,:,1) 
    500          trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:)      * zrfact2 * tmask(:,:,:) 
    501          trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:)      * zrfact2 * tmask(:,:,:) 
    502          trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:)      * zrfact2 * tmask(:,:,:) 
    503          trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:)      * zrfact2 * tmask(:,:,:) 
    504          trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d  (:,:,:)                * tmask(:,:,:) 
    505          trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3  (:,:,:)                * tmask(:,:,:) 
    506          trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4  (:,:,:)                * tmask(:,:,:) 
    507 # endif 
    508         ! 
    509598      ENDIF 
     599 
    510600      ! 
    511601      IF(ln_ctl)   THEN  ! print mean trends (used for debugging) 
     
    663753      END DO 
    664754      ! 
     755      ik100 = 10        !  last level where depth less than 100 m 
     756      DO jk = jpkm1, 1, -1 
     757         IF( gdept_1d(jk) > 100. )  iksed = jk - 1 
     758      END DO 
     759      IF (lwp) WRITE(numout,*) 
     760      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1 
     761      IF (lwp) WRITE(numout,*) 
     762      ! 
    665763      t_oce_co2_exp = 0._wp 
    666764      ! 
     
    702800      ztraz(:,:,:) = 0.e0 
    703801      zakz (:,:,:) = 0.e0 
    704       ztrb (:,:,:) = trn(:,:,:,jp_tra) 
     802      ztrb (:,:,:) = trb(:,:,:,jp_tra) 
    705803 
    706804      DO jk = 1, jpkm1 
     
    717815         !  first guess of the slopes interior values 
    718816         DO jk = 2, jpkm1 
    719             ztraz(:,:,jk) = ( trn(:,:,jk-1,jp_tra) - trn(:,:,jk,jp_tra) ) * tmask(:,:,jk) 
     817            ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk) 
    720818         END DO 
    721819         ztraz(:,:,1  ) = 0.0 
     
    748846                  zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1) 
    749847                  zew   = zwsink2(ji,jj,jk+1) 
    750                   psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 
     848                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 
    751849               END DO 
    752850            END DO 
     
    761859               DO ji = 1, jpi 
    762860                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 
    763                   trn(ji,jj,jk,jp_tra) = trn(ji,jj,jk,jp_tra) + zflx 
     861                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx 
    764862               END DO 
    765863            END DO 
     
    777875      END DO 
    778876 
    779       trn(:,:,:,jp_tra) = ztrb(:,:,:) 
     877      trb(:,:,:,jp_tra) = ztrb(:,:,:) 
    780878      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:) 
    781879      ! 
Note: See TracChangeset for help on using the changeset viewer.