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 10222 for NEMO/trunk/src/TOP/PISCES/SED/seddta.F90 – NEMO

Ignore:
Timestamp:
2018-10-25T11:42:23+02:00 (6 years ago)
Author:
aumont
Message:

update of the sediment module + enhancement, bug correction in PISCES generating O2 negative values

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/TOP/PISCES/SED/seddta.F90

    r7646 r10222  
    44   !! Sediment data  :  read sediment input data from a file 
    55   !!===================================================================== 
    6 #if defined key_sed 
     6 
    77   !! * Modules used 
    88   USE sed 
    99   USE sedarr 
    1010   USE iom 
     11   USE lib_mpp         ! distribued memory computing library 
    1112 
    1213   IMPLICIT NONE 
     
    1718 
    1819   !! *  Module variables 
    19    REAL(wp), DIMENSION(:), ALLOCATABLE ::  smask       ! mask for sediments points 
    2020   REAL(wp) ::  rsecday  ! number of second per a day 
    21    REAL(wp) ::  conv1    ! [m/day]--->[cm/s]   
    2221   REAL(wp) ::  conv2    ! [kg/m2/month]-->[g/cm2/s] ( 1 month has 30 days ) 
    23  
    24    INTEGER ::  numbio 
    25  
    26 #if defined key_sed_off 
    27    INTEGER ::  numoce 
    28 #endif 
    2922 
    3023   !! $Id$ 
     
    5447      INTEGER  ::  ji, jj, js, jw, ikt 
    5548 
    56       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zdta 
    57       REAL(wp), DIMENSION(:)  , ALLOCATABLE :: zdtap, zdtag 
    58  
     49      REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag 
     50      REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3, zwscal 
     51      REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep 
    5952 
    6053      !---------------------------------------------------------------------- 
     
    6457      !------------------------------------------------------------- 
    6558 
    66       WRITE(numsed,*) 
    67       WRITE(numsed,*) ' sed_dta : Bottom layer fields' 
    68       WRITE(numsed,*) ' ~~~~~~' 
    69       WRITE(numsed,*) ' Data from SMS model' 
    70       WRITE(numsed,*) 
     59      IF( ln_timing )  CALL timing_start('sed_dta') 
     60 
     61      IF (lwp) THEN 
     62         WRITE(numsed,*) 
     63         WRITE(numsed,*) ' sed_dta : Bottom layer fields' 
     64         WRITE(numsed,*) ' ~~~~~~' 
     65         WRITE(numsed,*) ' Data from SMS model' 
     66         WRITE(numsed,*) 
     67      ENDIF 
    7168 
    7269 
    7370      ! open file 
    7471      IF( kt == nitsed000 ) THEN 
    75          WRITE(numsed,*) ' sed_dta : Sediment fields' 
    76          CALL iom_open ( 'data_bio_bot'     , numbio ) 
    77 #if defined key_sed_off 
    78          CALL iom_open( 'data_oce_bot', numoce) 
    79 #endif 
     72         IF (lwp) WRITE(numsed,*) ' sed_dta : Sediment fields' 
     73         dtsed = r2dttrc 
    8074         rsecday = 60.* 60. * 24. 
    81          conv1   = 1.0e+2 / rsecday  
    82          conv2   = 1.0e+3 / ( 1.0e+4 * rsecday * 30. )  
    83  
    84          ! Compute sediment mask 
    85          ALLOCATE( zdta(jpi,jpj) )  
     75!         conv2   = 1.0e+3 / ( 1.0e+4 * rsecday * 30. ) 
     76         conv2 = 1.0e+3 /  1.0e+4  
     77         rdtsed(2:jpksed) = dtsed / ( denssol * por1(2:jpksed) ) 
     78      ENDIF 
     79 
     80      ! Initialization of temporaries arrays   
     81      zdtap(:)    = 0.  
     82      zdtag(:)    = 0.   
     83 
     84      ! reading variables 
     85      IF (lwp) WRITE(numsed,*) 
     86      IF (lwp) WRITE(numsed,*) ' sed_dta : Bottom layer fields at time  kt = ', kt 
     87      ! reading variables 
     88      ! 
     89      !    Sinking speeds of detritus is increased with depth as shown 
     90      !    by data and from the coagulation theory 
     91      !    ----------------------------------------------------------- 
     92      IF (ln_sediment_offline) THEN 
    8693         DO jj = 1, jpj 
    8794            DO ji = 1, jpi 
    88                ikt = MAX( INT( sbathy(ji,jj) ) - 1, 1 ) 
    89                zdta(ji,jj) = tmask(ji,jj,ikt)  
    90             ENDDO 
    91          ENDDO 
    92          ALLOCATE( smask(jpoce) ) 
    93          smask(:) = 0. 
    94          CALL pack_arr( jpoce, smask(1:jpoce), zdta(1:jpi,1:jpj), iarroce(1:jpoce) ) 
    95       ENDIF 
    96  
    97       ! Initialization of temporaries arrays   
    98       ALLOCATE( zdtap(jpoce) )    ;   zdtap(:)    = 0.  
    99       ALLOCATE( zdtag(jpoce) )    ;   zdtag(:)    = 0.   
    100  
    101       IF( MOD( kt - 1, nfreq ) == 0 ) THEN 
    102          ! reading variables 
    103          WRITE(numsed,*) 
    104          WRITE(numsed,*) ' sed_dta : Bottom layer fields at time  kt = ', kt 
    105          ! reading variables 
    106          trc_data(:,:,:) = 0. 
    107 #if ! defined key_sed_off 
    108          DO jj = 1,jpj 
     95               ikt = mbkt(ji,jj) 
     96               zwsbio4(ji,jj) = wsbio2 / rday 
     97               zwsbio3(ji,jj) = wsbio  / rday 
     98               zwscal (ji,jj) = wsbio2 / rday 
     99            END DO 
     100         END DO 
     101      ELSE 
     102         DO jj = 1, jpj 
    109103            DO ji = 1, jpi 
    110104               ikt = mbkt(ji,jj) 
    111                IF ( tmask(ji,jj,ikt) == 1 ) THEN 
    112                   trc_data(ji,jj,1)  = trn  (ji,jj,ikt,jptal) 
    113                   trc_data(ji,jj,2)  = trn  (ji,jj,ikt,jpdic) 
    114                   trc_data(ji,jj,3)  = trn  (ji,jj,ikt,jpno3) / 7.6 
    115                   trc_data(ji,jj,4)  = trn  (ji,jj,ikt,jppo4) / 122. 
    116                   trc_data(ji,jj,5)  = trn  (ji,jj,ikt,jpoxy) 
    117                   trc_data(ji,jj,6)  = trn  (ji,jj,ikt,jpsil) 
    118                   trc_data(ji,jj,7 ) = sinksil (ji,jj,ikt) 
    119                   trc_data(ji,jj,8 ) = sinking (ji,jj,ikt) 
    120                   trc_data(ji,jj,9 ) = sinking2(ji,jj,ikt) 
    121                   trc_data(ji,jj,10) = sinkcal (ji,jj,ikt) 
    122                   trc_data(ji,jj,11) = tsn     (ji,jj,ikt,jp_tem) 
    123                   trc_data(ji,jj,12) = tsn     (ji,jj,ikt,jp_sal) 
    124                ENDIF 
    125             ENDDO 
     105               zdep = e3t_n(ji,jj,ikt) / r2dttrc 
     106               zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) / rday ) 
     107               zwscal (ji,jj) = MIN( 0.99 * zdep, wscal (ji,jj,ikt) / rday ) 
     108               zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) / rday ) 
     109            END DO 
     110         END DO 
     111      ENDIF 
     112 
     113      trc_data(:,:,:) = 0. 
     114      DO jj = 1,jpj 
     115         DO ji = 1, jpi 
     116            ikt = mbkt(ji,jj) 
     117            IF ( tmask(ji,jj,ikt) == 1 ) THEN 
     118               trc_data(ji,jj,1)   = trb(ji,jj,ikt,jpsil) 
     119               trc_data(ji,jj,2)   = trb(ji,jj,ikt,jpoxy) 
     120               trc_data(ji,jj,3)   = trb(ji,jj,ikt,jpdic) 
     121               trc_data(ji,jj,4)   = trb(ji,jj,ikt,jpno3) / 7.625 
     122               trc_data(ji,jj,5)   = trb(ji,jj,ikt,jppo4) / 122. 
     123               trc_data(ji,jj,6)   = trb(ji,jj,ikt,jptal) 
     124               trc_data(ji,jj,7)   = trb(ji,jj,ikt,jpnh4) / 7.625 
     125               trc_data(ji,jj,8)   = 0.0 
     126               trc_data(ji,jj,9)   = 28.0E-3 
     127               trc_data(ji,jj,10)  = trb(ji,jj,ikt,jpfer) 
     128               trc_data(ji,jj,11 ) = MIN(trb(ji,jj,ikt,jpgsi), 1E-4) * zwsbio4(ji,jj) * 1E3 
     129               trc_data(ji,jj,12 ) = MIN(trb(ji,jj,ikt,jppoc), 1E-4) * zwsbio3(ji,jj) * 1E3 
     130               trc_data(ji,jj,13 ) = MIN(trb(ji,jj,ikt,jpgoc), 1E-4) * zwsbio4(ji,jj) * 1E3 
     131               trc_data(ji,jj,14)  = MIN(trb(ji,jj,ikt,jpcal), 1E-4) * zwscal (ji,jj) * 1E3 
     132               trc_data(ji,jj,15)  = tsn(ji,jj,ikt,jp_tem) 
     133               trc_data(ji,jj,16)  = tsn(ji,jj,ikt,jp_sal) 
     134               trc_data(ji,jj,17 ) = ( trb(ji,jj,ikt,jpsfe) * zwsbio3(ji,jj) + trb(ji,jj,ikt,jpbfe)  & 
     135               &                     * zwsbio4(ji,jj)  ) * 1E3 / ( trc_data(ji,jj,12 ) + trc_data(ji,jj,13 ) + rtrn ) 
     136               trc_data(ji,jj,17 ) = MIN(1E-3, trc_data(ji,jj,17 ) ) 
     137            ENDIF 
    126138         ENDDO 
    127  
    128 #else 
    129          CALL iom_get( numbio, jpdom_data, 'ALKBOT'     , trc_data(:,:,1 ) ) 
    130          CALL iom_get( numbio, jpdom_data, 'DICBOT'     , trc_data(:,:,2 ) ) 
    131          CALL iom_get( numbio, jpdom_data, 'NO3BOT'     , trc_data(:,:,3 ) ) 
    132          CALL iom_get( numbio, jpdom_data, 'PO4BOT'     , trc_data(:,:,4 ) ) 
    133          CALL iom_get( numbio, jpdom_data, 'O2BOT'      , trc_data(:,:,5 ) ) 
    134          CALL iom_get( numbio, jpdom_data, 'SIBOT'      , trc_data(:,:,6 ) ) 
    135          CALL iom_get( numbio, jpdom_data, 'OPALFLXBOT' , trc_data(:,:,7 ) )  
    136          CALL iom_get( numbio, jpdom_data, 'POCFLXBOT'  , trc_data(:,:,8 ) )  
    137          CALL iom_get( numbio, jpdom_data, 'GOCFLXBOT'  , trc_data(:,:,9 ) )  
    138          CALL iom_get( numbio, jpdom_data, 'CACO3FLXBOT', trc_data(:,:,10) )  
    139          CALL iom_get( numoce, jpdom_data, 'TBOT'       , trc_data(:,:,11) )  
    140          CALL iom_get( numoce, jpdom_data, 'SBOT'       , trc_data(:,:,12) )  
    141 #endif 
    142  
    143          ! Pore water initial concentration [mol/l] in  k=1 
    144          !------------------------------------------------- 
    145  
    146           ! Alkalinity ( 1 umol = 10-6equivalent ) 
    147          CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwalk), trc_data(1:jpi,1:jpj,1), iarroce(1:jpoce) ) 
    148          ! DIC  
    149          CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwdic), trc_data(1:jpi,1:jpj,2), iarroce(1:jpoce) ) 
    150          ! Nitrates (1 umol/l = 10-6 mol/l) 
    151          CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwno3), trc_data(1:jpi,1:jpj,3), iarroce(1:jpoce) ) 
    152          ! Phosphates (1 umol/l = 10-6 mol/l) 
    153          CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwpo4), trc_data(1:jpi,1:jpj,4), iarroce(1:jpoce) ) 
    154          ! Oxygen (1 umol/l = 10-6 mol/l) 
    155          CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwoxy), trc_data(1:jpi,1:jpj,5), iarroce(1:jpoce) )         
    156          ! Silicic Acid [mol.l-1] 
    157          CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwsil), trc_data(1:jpi,1:jpj,6), iarroce(1:jpoce) )                   
    158          ! DIC13 (mol/l)obtained from dc13 and DIC (12) and PDB  
    159          CALL iom_get ( numbio,jpdom_data,'DC13',zdta(:,:) ) 
    160          CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwc13), zdta(1:jpi,1:jpj), iarroce(1:jpoce) ) 
    161          pwcp_dta(1:jpoce,jwc13) = pdb * ( pwcp_dta(1:jpoce,jwc13) * 1.0e-3 + 1.0 )  & 
    162             &                          *   pwcp_dta(1:jpoce,jwdic)          
    163           
    164          !  Solid components :  
    165          !----------------------- 
    166          !  Sinking fluxes for OPAL in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
    167          CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,7), iarroce(1:jpoce) )  
    168          rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4 
    169          !  Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
    170          CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,8) , iarroce(1:jpoce) )       
    171          CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,9) , iarroce(1:jpoce) ) 
    172          rainrm_dta(1:jpoce,jspoc) =   ( zdtap(1:jpoce) +  zdtag(1:jpoce) ) * 1e-4 
    173          !  Sinking fluxes for Calcite in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
    174          CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,10), iarroce(1:jpoce) ) 
    175          rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4 
    176          ! vector temperature [°C] and salinity  
    177          CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,11), iarroce(1:jpoce) ) 
    178          CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,12), iarroce(1:jpoce) ) 
    179          
    180          ! Clay rain rate in [mol/(cm**2.s)]  
    181          ! inputs data in [kg.m-2.mois-1] ---> 1e+3/(1e+4*60*24*60*60) [g.cm-2.s-1]    
    182          ! divided after by molecular weight g.mol-1       
    183          zdta(:,:)   = 0. 
    184          CALL iom_get( numbio, jpdom_data, 'CLAY', zdta(:,:) )  
    185          CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsclay) , zdta(1:jpi,1:jpj), iarroce(1:jpoce) )       
    186          rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay) 
    187  
    188       ENDIF 
     139      ENDDO 
     140 
     141      ! Pore water initial concentration [mol/l] in  k=1 
     142      !------------------------------------------------- 
     143      DO jw = 1, jpwat 
     144         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jw), trc_data(1:jpi,1:jpj,jw), iarroce(1:jpoce) ) 
     145      END DO 
     146      !  Solid components :  
     147      !----------------------- 
     148      !  Sinking fluxes for OPAL in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
     149      CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,11), iarroce(1:jpoce) )  
     150      rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4 
     151      !  Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
     152      CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,12) , iarroce(1:jpoce) )       
     153      CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,13) , iarroce(1:jpoce) ) 
     154      DO ji = 1, jpoce 
     155!        zkapp  = MIN( (1.0 - 0.02 ) * reac_poc, 3731.0 * max(100.0, zkbot(ji) )**(-1.011) / ( 365.0 * 24.0 * 3600.0 ) ) 
     156!        zkapp   = MIN( 0.98 * reac_poc, 100.0 * max(100.0, zkbot(ji) )**(-0.6) / ( 365.0 * 24.0 * 3600.0 ) ) 
     157!        zratio = ( ( 1.0 - 0.02 ) * reac_poc + 0.02 * reac_poc * 0. - zkapp) / ( ( 0.02 - 1.0 ) * reac_poc / 100. - 0.02 * reac_poc * 0. + zkapp ) 
     158!        zf1    = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc ) 
     159!        zf1    = MIN(0.98, MAX(0., zf1 ) ) 
     160         zf1    = 0.48 
     161         zf0    = 1.0 - 0.02 - zf1 
     162         zf2    = 0.02 
     163         rainrm_dta(ji,jspoc) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf0 
     164         rainrm_dta(ji,jspos) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf1 
     165         rainrm_dta(ji,jspor) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf2 
     166      END DO 
     167      !  Sinking fluxes for Calcite in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
     168      CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,14), iarroce(1:jpoce) ) 
     169      rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4 
     170      ! vector temperature [°C] and salinity  
     171      CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) ) 
     172      CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) ) 
     173       
     174      ! Clay rain rate in [mol/(cm**2.s)]  
     175      ! inputs data in [kg.m-2.sec-1] ---> 1e+3/(1e+4) [g.cm-2.s-1]    
     176      ! divided after by molecular weight g.mol-1       
     177      CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jsclay), dust(1:jpi,1:jpj), iarroce(1:jpoce) ) 
     178      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay)   & 
     179      &                            + wacc(1:jpoce) * por1(2) * denssol / mol_wgt(jsclay) / ( rsecday * 365.0 ) 
     180      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * 0.965 
     181      rainrm_dta(1:jpoce,jsfeo)  = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 / 0.965 
     182!    rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay) 
     183 
     184      ! Iron monosulphide rain rates. Set to 0 
     185      rainrm_dta(1:jpoce,jsfes)  = 0.  
     186 
     187      ! Fe/C ratio in sinking particles that fall to the sediments 
     188      CALL pack_arr ( jpoce,  fecratio(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) ) 
     189 
     190      sedligand(:,1) = 1.E-9 
    189191 
    190192      ! sediment pore water at 1st layer (k=1) 
    191193      DO jw = 1, jpwat 
    192          pwcp(1:jpoce,1,jw) = pwcp_dta(1:jpoce,jw) * smask(1:jpoce) 
     194         pwcp(1:jpoce,1,jw) = pwcp_dta(1:jpoce,jw) 
    193195      ENDDO 
    194196 
    195197      !  rain 
    196198      DO js = 1, jpsol 
    197          rainrm(1:jpoce,js) = rainrm_dta(1:jpoce,js) * smask(1:jpoce) 
     199         rainrm(1:jpoce,js) = rainrm_dta(1:jpoce,js) 
    198200      ENDDO 
    199201 
     
    212214      dzdep(1:jpoce) = raintg(1:jpoce) * rdtsed(2)  
    213215 
    214  
    215       DEALLOCATE( zdta )  
    216       DEALLOCATE( zdtap    ) ;  DEALLOCATE( zdtag    )  
    217  
    218       IF( kt == nitsedend )   THEN 
    219          CALL iom_close ( numbio ) 
    220 #if defined key_sed_off 
    221          CALL iom_close ( numoce ) 
    222 #endif 
    223       ENDIF 
     216      IF( lk_iomput ) THEN 
     217          IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * conv2 * 1E4 ) 
     218          IF( iom_use("sflxcal" ) )  CALL iom_put( "sflxcal", trc_data(:,:,13) ) 
     219          IF( iom_use("sflxbsi" ) )  CALL iom_put( "sflxbsi", trc_data(:,:,10) ) 
     220          IF( iom_use("sflxpoc" ) )  CALL iom_put( "sflxpoc", trc_data(:,:,11) + trc_data(:,:,12) ) 
     221      ENDIF 
     222 
     223      IF( ln_timing )  CALL timing_stop('sed_dta') 
    224224       
    225225   END SUBROUTINE sed_dta 
    226226 
    227 #else 
    228    !!====================================================================== 
    229    !! MODULE seddta  :   Dummy module  
    230    !!====================================================================== 
    231    !! $Id$ 
    232 CONTAINS 
    233    SUBROUTINE sed_dta ( kt ) 
    234      INTEGER, INTENT(in) :: kt 
    235      WRITE(*,*) 'sed_stp: You should not have seen this print! error?', kt  
    236   END SUBROUTINE sed_dta 
    237 #endif 
    238  
    239227END MODULE seddta 
Note: See TracChangeset for help on using the changeset viewer.