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 6291 for trunk – NEMO

Changeset 6291 for trunk


Ignore:
Timestamp:
2016-02-05T05:51:35+01:00 (8 years ago)
Author:
cetlod
Message:

trunk : bugfix on PISCES carbonate chemistry, see ticket #1672

Location:
trunk/NEMOGCM/NEMO/TOP_SRC/PISCES
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90

    r6140 r6291  
    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 
    34    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2    ! Solubilities of O2 and CO2 
     33   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 
    3535 
    3636   REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm 
     
    7676   REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate 
    7777   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  
    9178 
    9279   REAL(wp) ::   ft1    =    0.000067   ! constants for calculate concentrations for fluorides 
    9380   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 
    12081 
    12182   !                                    ! volumetric solubility constants for o2 in ml/L   
     
    196157         DO ji = 1, jpi 
    197158            !                             ! SET ABSOLUTE TEMPERATURE 
    198             ztkel = tsn(ji,jj,1,jp_tem) + 273.16 
     159            ztkel = tsn(ji,jj,1,jp_tem) + 273.15 
    199160            zt    = ztkel * 0.01 
    200161            zt2   = zt * zt 
     
    205166            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
    206167            zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 ) 
    207             !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 in ml/L (EQ. 8, GARCIA AND GORDON, 1992) 
    208             ztgg  = LOG( ( 298.15 - tsn(ji,jj,1,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
    209             ztgg2 = ztgg  * ztgg 
    210             ztgg3 = ztgg2 * ztgg 
    211             ztgg4 = ztgg3 * ztgg 
    212             ztgg5 = ztgg4 * ztgg 
    213             zoxy  = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5   & 
    214                    + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) +  ox10 * zsal2 
    215  
    216168            !                             ! SET SOLUBILITIES OF O2 AND CO2  
    217             chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm) 
    218             chemc(ji,jj,2) = ( EXP( zoxy  ) * o2atm ) * oxyco              ! mol/(L atm) 
     169            chemc(ji,jj) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm) 
    219170            ! 
    220171         END DO 
     
    226177         DO jj = 1, jpj 
    227178            DO ji = 1, jpi 
    228               ztkel = tsn(ji,jj,jk,jp_tem) + 273.16 
     179              ztkel = tsn(ji,jj,jk,jp_tem) + 273.15 
    229180              zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35. 
    230181              zsal2 = zsal * zsal 
     
    242193 
    243194 
     195 
    244196      ! CHEMICAL CONSTANTS - DEEP OCEAN 
    245197      ! ------------------------------- 
     
    252204 
    253205               ! SET ABSOLUTE TEMPERATURE 
    254                ztkel   = tsn(ji,jj,jk,jp_tem) + 273.16 
     206               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.15 
    255207               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35. 
    256208               zsqrt  = SQRT( zsal ) 
     
    273225 
    274226               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 
    275                zcks    = EXP(  ks1 * ztr + ks0 + ks2 * zlogt                           & 
    276                   &                     + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt   & 
    277                   &                     + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis      & 
    278                   &                     + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal )  ) 
     227               zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         & 
     228               &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt & 
     229               &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    & 
     230               &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         & 
     231               &         + LOG(1.0 - 0.001005 * zsal)) 
     232               ! 
     233               aphscale(ji,jj,jk) = ( 1. + zst / zcks ) 
    279234 
    280235               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 
    281                zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  ) 
     236               zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   & 
     237               &         + LOG(1.0d0 - 0.001005d0*zsal)            & 
     238               &         + LOG(1.0d0 + zst/zcks)) 
    282239 
    283240               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 
    284                zckb    = ( cb0 + cb1 * zsqrt + cb2  * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr   & 
    285                   &    + ( cb5 + cb6 * zsqrt + cb7  * zsal )                                            & 
    286                   &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             & 
    287                   &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  ) 
     241               zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        & 
     242               &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         & 
     243               &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   & 
     244               &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      &  
     245               &      * zlogt + 0.053105*zsqrt*ztkel 
     246 
    288247 
    289248               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal 
     
    291250 
    292251               ! PKW (H2O) (DICKSON AND RILEY, 1979) 
    293                zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal 
    294  
     252               zckw = -13847.26*ztr + 148.9652 - 23.6521 * zlogt    &  
     253               &     + (118.67*ztr - 5.977 + 1.0495 * zlogt)        & 
     254               &     * zsqrt - 0.01615 * zsal 
    295255 
    296256               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
     
    367327      !!                     ***  ROUTINE p4z_che_alloc  *** 
    368328      !!---------------------------------------------------------------------- 
    369       ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,2), chemo2(jpi,jpj,jpk), STAT=p4z_che_alloc ) 
     329      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj), chemo2(jpi,jpj,jpk),   & 
     330      &         STAT=p4z_che_alloc ) 
    370331      ! 
    371332      IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90

    r6140 r6291  
    8282      ! 
    8383      INTEGER  ::   ji, jj, jm, iind, iindm1 
    84       REAL(wp) ::   ztc, ztc2, ztc3, zws, zkgwan 
     84      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan 
    8585      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact 
    8686      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 
     
    131131 
    132132               ! CALCULATE [ALK]([CO3--], [HCO3-]) 
    133                zalk  = zalka - (  akw3(ji,jj,1) / zph - zph + zbot / ( 1.+ zph / akb3(ji,jj,1) )  ) 
     133               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph / aphscale(ji,jj,1)    & 
     134               &       + zbot / ( 1.+ zph / akb3(ji,jj,1) )  ) 
    134135 
    135136               ! CALCULATE [H+] AND [H2CO3] 
     
    156157            ztc2 = ztc * ztc 
    157158            ztc3 = ztc * ztc2  
     159            ztc4 = ztc2 * ztc2  
    158160            ! Compute the schmidt Number both O2 and CO2 
    159             zsch_co2 = 2073.1 - 125.62 * ztc + 3.6276 * ztc2 - 0.043126 * ztc3 
    160             zsch_o2  = 1953.4 - 128.0  * ztc + 3.9918 * ztc2 - 0.050091 * ztc3 
     161            zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4 
     162            zsch_o2  = 1920.4 - 135.6  * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4 
    161163            !  wind speed  
    162164            zws  = wndm(ji,jj) * wndm(ji,jj) 
    163165            ! Compute the piston velocity for O2 and CO2 
    164             zkgwan = 0.3 * zws  + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946  * ztc2 ) 
     166            zkgwan = 0.251 * zws 
    165167            zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 
    166168# if defined key_degrad 
     
    176178         DO ji = 1, jpi 
    177179            ! Compute CO2 flux for the sea and air 
    178             zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)   ! (mol/L) * (m/s) 
     180            zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj) * zkgco2(ji,jj)   ! (mol/L) * (m/s) 
    179181            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ? 
    180182            oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 
     
    183185 
    184186            ! Compute O2 flux  
    185             zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
     187            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
    186188            zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj) 
    187189            zoflx(ji,jj) = zfld16 - zflu16 
     
    216218         ENDIF 
    217219         IF( iom_use( "Dpco2" ) ) THEN 
    218            zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
     220           zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1) 
    219221           CALL iom_put( "Dpco2" ,  zw2d ) 
    220222         ENDIF 
    221223         IF( iom_use( "Dpo2" ) )  THEN 
    222            zw2d(:,:) = ( atcox * patm(:,:) - trb(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1) 
     224           zw2d(:,:) = ( atcox * patm(:,:) - atcox * trn(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) 
    223225           CALL iom_put( "Dpo2"  , zw2d ) 
    224226         ENDIF 
     
    232234            trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1)  
    233235            trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1)  
    234             trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)  
     236            trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1)  
    235237         ENDIF 
    236238      ENDIF 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

    r5836 r6291  
    8888                  zalka = trb(ji,jj,jk,jptal) / zfact 
    8989                  ! CALCULATE [ALK]([CO3--], [HCO3-]) 
    90                   zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 
     90                  zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph / ( aphscale(ji,jj,jk) + rtrn )  & 
     91                  &       + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 
    9192                  ! CALCULATE [H+] and [CO3--] 
    9293                  zaldi = zdic - zalk 
     
    149150         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) ) 
    150151      ELSE 
    151          trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:) 
    152          trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:) 
    153          trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:) 
     152         IF( ln_diatrc ) THEN 
     153            trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:) 
     154            trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:) 
     155            trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:) 
     156         ENDIF 
    154157      ENDIF 
    155158      ! 
  • trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90

    r5385 r6291  
    101101   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi         !: ??? 
    102102   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   excess     !: ??? 
     103   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aphscale   !:  
     104 
    103105 
    104106   !!* Temperature dependancy of SMS terms 
     
    154156         &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,       & 
    155157         &      akw3(jpi,jpj,jpk)    , borat (jpi,jpj,jpk) ,       & 
    156          &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,     STAT=ierr(4) ) 
     158         &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,       & 
     159         &      aphscale(jpi,jpj,jpk),                           STAT=ierr(4) ) 
    157160         ! 
    158161      !* Temperature dependancy of SMS terms 
Note: See TracChangeset for help on using the changeset viewer.