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 7607 for branches/2015/nemo_v3_6_STABLE – NEMO

Ignore:
Timestamp:
2017-01-25T16:37:31+01:00 (7 years ago)
Author:
cetlod
Message:

v3.6 stable : add missing features for CMIP6 exercise, see ticket #1834

Location:
branches/2015/nemo_v3_6_STABLE/NEMOGCM
Files:
21 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_cfg

    r4690 r7607  
    11!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    2 !! NEMO/LIM-2 : Ice configuration namelist. Overwrites SHARED/namelist_ice_lim2_ref 
     2!! NEMO/LIM-3 : Ice configuration namelist. Overwrites SHARED/namelist_ice_lim3_ref 
    33!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    44 
     
    1515!----------------------------------------------------------------------- 
    1616/ 
     17!------------------------------------------------------------------------------ 
     18&namicehdf     !   Ice horizontal diffusion 
     19!------------------------------------------------------------------------------ 
     20   nn_ahi0        =  2          !  horizontal diffusivity computation 
     21/ 
    1722!----------------------------------------------------------------------- 
    1823&namicethd     !   ice thermodynamic 
     
    2429/ 
    2530!----------------------------------------------------------------------- 
    26 &namiceitdme   !   parameters for mechanical redistribution of ice  
     31&namiceitdme   !   parameters for mechanical redistribution of ice 
    2732!----------------------------------------------------------------------- 
    2833/ 
     
    3136!----------------------------------------------------------------------- 
    3237/ 
     38!------------------------------------------------------------------------------ 
     39&namiceitd     !   Ice discretization 
     40!------------------------------------------------------------------------------ 
     41/ 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/field_def.xml

    r7511 r7607  
    225225         <field id="icealb_cea"   long_name="Ice albedo (cell average)"                 standard_name="sea_ice_albedo"                            unit="1"        /> 
    226226         <field id="calving_cea"  long_name="Calving"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s"  /> 
     227         <field id="iceberg_cea"  long_name="Iceberg"                                   standard_name="water_flux_into_sea_water_from_icebergs"   unit="kg/m2/s"  /> 
     228         <field id="iceshelf_cea"  long_name="Iceshelf"                                   standard_name="water_flux_into_sea_water_from_iceshelf"   unit="kg/m2/s"  /> 
    227229 
    228230         <!-- available if key_oasis3 + conservative method --> 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r6316 r7607  
    6969   rn_relast      =    0.333       !  ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast  
    7070                                   !     advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 
    71    nn_ahi0        =    2           !  horizontal diffusivity computation 
    72                                    !     0: use rn_ahi0_ref 
    73                                    !     1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 
    74                                    !     2: use rn_ahi0_ref x grid cell length      / ( 2deg mean grid cell length ) 
    75    rn_ahi0_ref    = 350.0          !  horizontal sea ice diffusivity (m2/s)  
    76                                    !     if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 
    7771/ 
    7872!------------------------------------------------------------------------------ 
    7973&namicehdf     !   Ice horizontal diffusion 
    8074!------------------------------------------------------------------------------ 
     75   nn_ahi0        =    -1          !  horizontal diffusivity computation 
     76                                   !    -1: no diffusion (bypass limhdf) 
     77                                   !     0: use rn_ahi0_ref 
     78                                   !     1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 
     79                                   !     2: use rn_ahi0_ref x grid cell length      / ( 2deg mean grid cell length ) 
     80   rn_ahi0_ref    = 350.0          !  horizontal sea ice diffusivity (m2/s) 
     81                                   !     if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 
    8182   nn_convfrq     = 5              !  convergence check frequency of the Crant-Nicholson scheme (perf. optimization) 
    8283/ 
     
    9899                                   !     0: k = k0 + beta.S/T (Untersteiner, 1964) 
    99100                                   !     1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) 
     101   rn_cdsn     = 0.31              !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 
     102                                   !  Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    100103   nn_monocat  = 0                 !  virtual ITD mono-category parameterizations (1, jpl = 1 only) or not (0) 
    101104                                   !     2: simple piling instead of ridging --- temporary option 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/SHARED/namelist_ref

    r6319 r7607  
    500500&namsbc_alb    !   albedo parameters 
    501501!----------------------------------------------------------------------- 
    502    nn_ice_alb  =    0   !  parameterization of ice/snow albedo 
    503                         !     0: Shine & Henderson-Sellers (JGR 1985) 
    504                         !     1: "home made" based on Brandt et al. (J. Climate 2005) 
    505                         !                         and Grenfell & Perovich (JGR 2004) 
    506    rn_albice   =  0.53  !  albedo of bare puddled ice (values from 0.49 to 0.58) 
    507                         !     0.53 (default) => if nn_ice_alb=0 
    508                         !     0.50 (default) => if nn_ice_alb=1 
     502   nn_ice_alb   =    1   !  parameterization of ice/snow albedo 
     503                         !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo 
     504                         !     1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004), 
     505                         !        giving cloud-sky albedo 
     506   rn_alb_sdry  =  0.85  !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky) 
     507   rn_alb_smlt  =  0.75  !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' ) 
     508   rn_alb_idry  =  0.60  !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' ) 
     509   rn_alb_imlt  =  0.50  !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' ) 
    509510/ 
    510511!----------------------------------------------------------------------- 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r6963 r7607  
    212212   REAL(wp), PUBLIC ::   rn_betas         !: coef. for partitioning of snowfall between leads and sea ice 
    213213   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
     214   REAL(wp), PUBLIC ::   rn_cdsn          !: thermal conductivity of the snow [W/m/K] 
    214215   REAL(wp), PUBLIC ::   nn_conv_dif      !: maximal number of iterations for heat diffusion 
    215216   REAL(wp), PUBLIC ::   rn_terr_dif      !: maximal tolerated error (C) for heat diffusion 
     
    320321   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i     !: Sea-Ice Age (days) 
    321322   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i    !: Sea-Ice Age times ice area (days) 
     323 
    322324   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i    !: brine volume 
    323325 
     
    406408   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s]  
    407409   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s]  
     410 
    408411   ! 
    409412   !!---------------------------------------------------------------------- 
     
    463466         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     & 
    464467         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     & 
    465          &      om_i (jpi,jpj)                              , STAT=ierr(ii) ) 
     468         &      om_i (jpi,jpj) , STAT=ierr(ii) ) 
    466469      ii = ii + 1 
    467470      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r5123 r7607  
    244244      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    245245      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, & 
    246          &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 
    247       INTEGER  ::   ji, jj 
    248       REAL(wp) ::   za00, zd_max 
     246         &                nn_nevp, rn_relast 
    249247      !!------------------------------------------------------------------- 
    250248 
     
    272270         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp 
    273271         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast 
    274          WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0 
    275          WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref 
    276272      ENDIF 
    277273      ! 
     
    279275      rhoco  = rau0  * rn_cio 
    280276      ! 
    281       !  Diffusion coefficients 
    282       SELECT CASE( nn_ahi0 ) 
    283  
    284       CASE( 0 ) 
    285          ahiu(:,:) = rn_ahi0_ref 
    286          ahiv(:,:) = rn_ahi0_ref 
    287  
    288          IF(lwp) WRITE(numout,*) '' 
    289          IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref' 
    290  
    291       CASE( 1 )  
    292  
    293          zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
    294          IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain 
    295           
    296          ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
    297                                                         !                    (60° = min latitude for ice cover)   
    298          ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 
    299  
    300          IF(lwp) WRITE(numout,*) '' 
    301          IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 
    302          IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp  
    303           
    304       CASE( 2 )  
    305  
    306          zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
    307          IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
    308           
    309          za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
    310                                                  !                    (60° = min latitude for ice cover)   
    311          DO jj = 1, jpj 
    312             DO ji = 1, jpi 
    313                ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 
    314                ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 
    315             END DO 
    316          END DO 
    317          ! 
    318          IF(lwp) WRITE(numout,*) '' 
    319          IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1' 
    320          IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 
    321           
    322       END SELECT 
    323  
    324277   END SUBROUTINE lim_dyn_init 
    325278 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r6476 r7607  
    233233      !!------------------------------------------------------------------- 
    234234      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    235       NAMELIST/namicehdf/ nn_convfrq  
    236       !!------------------------------------------------------------------- 
    237       ! 
    238       IF(lwp) THEN 
    239          WRITE(numout,*) 
    240          WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 
    241          WRITE(numout,*) '~~~~~~~' 
    242       ENDIF 
     235      NAMELIST/namicehdf/  nn_ahi0, rn_ahi0_ref, nn_convfrq  
     236      INTEGER  ::   ji, jj 
     237      REAL(wp) ::   za00, zd_max 
     238      !!------------------------------------------------------------------- 
    243239      ! 
    244240      REWIND( numnam_ice_ref )              ! Namelist namicehdf in reference namelist : Ice horizontal diffusion 
     
    253249      IF(lwp) THEN                          ! control print 
    254250         WRITE(numout,*) 
    255          WRITE(numout,*)'   Namelist of ice parameters for ice horizontal diffusion computation ' 
    256          WRITE(numout,*)'      convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq 
     251         WRITE(numout,*) 'lim_hdf_init : Ice horizontal diffusion' 
     252         WRITE(numout,*) '~~~~~~~~~~~' 
     253         WRITE(numout,*) '   horizontal diffusivity calculation                          nn_ahi0      = ', nn_ahi0 
     254         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)                  rn_ahi0_ref  = ', rn_ahi0_ref 
     255         WRITE(numout,*) '   convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq 
    257256      ENDIF 
     257      ! 
     258      !  Diffusion coefficients 
     259      SELECT CASE( nn_ahi0 ) 
     260 
     261      CASE( -1 ) 
     262         ahiu(:,:) = 0._wp 
     263         ahiv(:,:) = 0._wp 
     264 
     265         IF(lwp) WRITE(numout,*) '' 
     266         IF(lwp) WRITE(numout,*) '   No sea-ice diffusion applied' 
     267 
     268      CASE( 0 ) 
     269         ahiu(:,:) = rn_ahi0_ref 
     270         ahiv(:,:) = rn_ahi0_ref 
     271 
     272         IF(lwp) WRITE(numout,*) '' 
     273         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref' 
     274 
     275      CASE( 1 )  
     276 
     277         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     278         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain 
     279          
     280         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 
     281                                                        !                    (60deg = min latitude for ice cover)   
     282         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 
     283 
     284         IF(lwp) WRITE(numout,*) '' 
     285         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 
     286         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp  
     287          
     288      CASE( 2 )  
     289 
     290         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     291         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
     292          
     293         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 
     294                                                 !                    (60deg = min latitude for ice cover)   
     295         DO jj = 1, jpj 
     296            DO ji = 1, jpi 
     297               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 
     298               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 
     299            END DO 
     300         END DO 
     301         ! 
     302         IF(lwp) WRITE(numout,*) '' 
     303         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1' 
     304         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 
     305          
     306      END SELECT 
    258307      ! 
    259308   END SUBROUTINE lim_hdf_init 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r6399 r7607  
    613613         CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
    614614         CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj ) 
    615          !          
    616615      END SELECT 
    617616 
     
    633632      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    634633      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       & 
    635          &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, & 
    636          &                nn_monocat, ln_it_qnsice 
     634         &                rn_himin, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon,         & 
     635         &                rn_cdsn, nn_monocat, ln_it_qnsice 
    637636      !!------------------------------------------------------------------- 
    638637      ! 
     
    673672         WRITE(numout,*)'      maximal err. on T for heat diffusion computation        rn_terr_dif  = ', rn_terr_dif 
    674673         WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon 
     674         WRITE(numout,*)'      thermal conductivity of the snow                        rn_cdsn      = ', rn_cdsn 
    675675         WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i 
    676676         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5512 r7607  
    376376 
    377377               ! Effective thickness he (zhe) 
    378                zfac     = 1._wp / ( rcdsn + zkimean ) 
    379                zratio_s = rcdsn   * zfac 
     378               zfac     = 1._wp / ( rn_cdsn + zkimean ) 
     379               zratio_s = rn_cdsn   * zfac 
    380380               zratio_i = zkimean * zfac 
    381381               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 
     
    400400         DO ji = kideb, kiut 
    401401            zfac                  =  1. / MAX( epsi10 , zh_s(ji) ) 
    402             zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac 
    403             zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac 
     402            zkappa_s(ji,0)        = zghe(ji) * rn_cdsn * zfac 
     403            zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cdsn * zfac 
    404404         END DO 
    405405 
    406406         DO jk = 1, nlay_s-1 
    407407            DO ji = kideb , kiut 
    408                zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 
     408               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 
    409409            END DO 
    410410         END DO 
     
    422422            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac 
    423423            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 
    424             zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / &  
    425            &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 
     424            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / &  
     425           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) ) 
    426426            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    427427         END DO 
     
    697697               &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
    698698         END DO 
     699 
    699700         ! 
    700701         !-------------------------------------------------------------------------- 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r5147 r7607  
    6565#if defined key_lim3 || defined key_cice 
    6666   REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3] 
    67    REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice 
    68    REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow 
    69    REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice  
     67   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K] 
     68   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat of fresh ice                            [J/kg/K] 
    7069   REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
    7170   REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
     
    8382   REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3] 
    8483   REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3] 
     84#endif 
     85#if defined key_cice 
     86   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K], now namelist parameter for LIM3 
    8587#endif 
    8688#if defined key_lim3 
     
    177179      IF(lwp) THEN 
    178180         WRITE(numout,*) 
     181#if defined key_cice 
    179182         WRITE(numout,*) '          thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K' 
    180          WRITE(numout,*) '          thermal conductivity of the ice           = ', rcdic   , ' J/s/m/K' 
     183#endif 
     184         WRITE(numout,*) '          thermal conductivity of pure ice          = ', rcdic   , ' J/s/m/K' 
    181185         WRITE(numout,*) '          fresh ice specific heat                   = ', cpic    , ' J/kg/K' 
    182186         WRITE(numout,*) '          latent heat of fusion of fresh ice / snow = ', lfus    , ' J/kg' 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

    r6323 r7607  
    3939   !                             !!* namelist namsbc_alb 
    4040   INTEGER  ::   nn_ice_alb 
    41    REAL(wp) ::   rn_albice 
     41   REAL(wp) ::   rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt 
    4242 
    4343   !!---------------------------------------------------------------------- 
     
    101101      IF( albd_init == 0 )   CALL albedo_init      ! initialization  
    102102 
     103      ralb_sf = rn_alb_sdry ! dry snow 
     104      ralb_sm = rn_alb_smlt ! melting snow 
     105      ralb_if = rn_alb_idry ! bare frozen ice 
     106      ralb_im = rn_alb_imlt ! bare puddled ice  
    103107       
    104108      SELECT CASE ( nn_ice_alb ) 
     
    109113      CASE( 0 ) 
    110114        
    111          ralb_sf = 0.80       ! dry snow 
    112          ralb_sm = 0.65       ! melting snow 
    113          ralb_if = 0.72       ! bare frozen ice 
    114          ralb_im = rn_albice  ! bare puddled ice  
    115           
     115         !ralb_sf = 0.80       ! dry snow 
     116         !ralb_sm = 0.65       ! melting snow 
     117         !ralb_if = 0.72       ! bare frozen ice 
     118         !ralb_im = ...        ! bare puddled ice  
     119 
    116120         !  Computation of ice albedo (free of snow) 
    117121         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im 
     
    163167      CASE( 1 )  
    164168 
    165          ralb_im = rn_albice  ! bare puddled ice 
     169!        ralb_im = ...        ! bare puddled ice 
    166170! compilation of values from literature 
    167          ralb_sf = 0.85      ! dry snow 
    168          ralb_sm = 0.75      ! melting snow 
    169          ralb_if = 0.60      ! bare frozen ice 
     171!        ralb_sf = 0.85      ! dry snow 
     172!        ralb_sm = 0.75      ! melting snow 
     173!        ralb_if = 0.60      ! bare frozen ice 
    170174! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 
    171175!         ralb_sf = 0.85       ! dry snow 
     
    248252      !!---------------------------------------------------------------------- 
    249253      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    250       NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice  
     254      NAMELIST/namsbc_alb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry , rn_alb_imlt 
    251255      !!---------------------------------------------------------------------- 
    252256      ! 
     
    268272         WRITE(numout,*) '   Namelist namsbc_alb : albedo ' 
    269273         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb 
    270          WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice 
     274         WRITE(numout,*) '      albedo of dry snow                                  rn_alb_sdry = ', rn_alb_sdry 
     275         WRITE(numout,*) '      albedo of melting snow                              rn_alb_smlt = ', rn_alb_smlt 
     276         WRITE(numout,*) '      albedo of dry ice                                   rn_alb_idry = ', rn_alb_idry 
     277         WRITE(numout,*) '      albedo of bare puddled ice                          rn_alb_imlt = ', rn_alb_imlt 
    271278      ENDIF 
    272279      ! 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r5407 r7607  
    113113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff        [Kg/m2/s]   
    114114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwfisf , fwfisf_b !: ice shelf melting   [Kg/m2/s]   
     115   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwficb , fwficb_b !: iceberg melting [Kg/m2/s]   
    115116   !! 
    116117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] jpi,jpj,jpts 
     
    164165         ! 
    165166      ALLOCATE( fwfisf  (jpi,jpj), rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     & 
    166          &      fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 
     167         &      fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) ,     & 
     168         &      fwficb  (jpi,jpj), fwficb_b(jpi,jpj), STAT=ierr(3) ) 
    167169         ! 
    168170      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r7494 r7607  
    4343   USE eosbn2 
    4444   USE sbcrnf   , ONLY : l_rnfcpl 
     45   USE sbcisf   , ONLY : l_isfcpl 
    4546#if defined key_cpl_carbon_cycle 
    4647   USE p4zflx, ONLY : oce_co2 
     
    105106   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness  
    106107   INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level 
    107    INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received 
     108   INTEGER, PARAMETER ::   jpr_isf    = 43 
     109   INTEGER, PARAMETER ::   jpr_icb    = 44 
     110   INTEGER, PARAMETER ::   jprcv      = 44            ! total number of fields received 
    108111 
    109112   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere 
     
    149152   ! Received from the atmosphere                     ! 
    150153   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 
    151    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                         
     154   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_icb, sn_rcv_isf                                
    152155   ! Other namelist parameters                        ! 
    153156   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
     
    219222         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
    220223         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   & 
    221          &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask 
     224         &                  sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf, nn_cplmodel  , ln_usecplmask 
    222225      !!--------------------------------------------------------------------- 
    223226      ! 
     
    258261         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')' 
    259262         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')' 
     263         WRITE(numout,*)'      iceberg                         = ', TRIM(sn_rcv_icb%cldes   ), ' (', TRIM(sn_rcv_icb%clcat   ), ')' 
     264         WRITE(numout,*)'      ice shelf                       = ', TRIM(sn_rcv_isf%cldes   ), ' (', TRIM(sn_rcv_isf%clcat   ), ')' 
    260265         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    261266         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')' 
     
    397402      END SELECT 
    398403 
    399       !                                                      ! ------------------------- ! 
    400       !                                                      !     Runoffs & Calving     !    
    401       !                                                      ! ------------------------- ! 
     404 
     405      !                                                      ! ---------------------------------------------------- ! 
     406      !                                                      !     Runoffs, Calving, Iceberg, Iceshelf cavities     !    
     407      !                                                      ! ---------------------------------------------------- ! 
    402408      srcv(jpr_rnf   )%clname = 'O_Runoff' 
    403409      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN 
     
    409415      ENDIF 
    410416      ! 
    411       srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE. 
     417      srcv(jpr_cal)%clname = 'OCalving'   ;  IF( TRIM( sn_rcv_cal%cldes) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE. 
     418      srcv(jpr_isf)%clname = 'OIcshelf'   ;  IF( TRIM( sn_rcv_isf%cldes) == 'coupled' )   srcv(jpr_isf)%laction = .TRUE. 
     419      srcv(jpr_icb)%clname = 'OIceberg'   ;  IF( TRIM( sn_rcv_icb%cldes) == 'coupled' )   srcv(jpr_icb)%laction = .TRUE. 
     420 
     421      IF( srcv(jpr_isf)%laction .AND. nn_isf > 0 ) THEN 
     422         l_isfcpl             = .TRUE.                      ! -> no need to read isf in sbcisf 
     423         IF(lwp) WRITE(numout,*) 
     424         IF(lwp) WRITE(numout,*) '   iceshelf received from oasis ' 
     425      ENDIF 
    412426 
    413427      !                                                      ! ------------------------- ! 
     
    10711085         ENDIF 
    10721086         ! 
     1087         !    
    10731088         !                                                        ! runoffs and calving (added in emp) 
    1074          IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1) 
     1089         IF( srcv(jpr_rnf)%laction )     rnf(:,:)  = frcv(jpr_rnf)%z3(:,:,1) 
    10751090         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1) 
     1091 
     1092         IF( srcv(jpr_icb)%laction )  THEN  
     1093             fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1) 
     1094             rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runfofs 
     1095         ENDIF 
     1096         IF( srcv(jpr_isf)%laction )  fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)   
    10761097          
    10771098         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:) 
     
    10911112            ENDIF 
    10921113         ENDIF 
     1114         ! 
     1115         IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting 
     1116         ! 
    10931117         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:) 
    10941118         ELSE                   ;   qns(:,:) =                              zqns(:,:) 
     
    14681492      ENDIF 
    14691493 
     1494      IF( srcv(jpr_icb)%laction )  THEN  
     1495         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1) 
     1496         rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runoffs 
     1497         CALL iom_put( 'iceberg_cea', frcv(jpr_icb)%z3(:,:,1) ) 
     1498      ENDIF 
     1499      IF( srcv(jpr_isf)%laction )  THEN 
     1500        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)   
     1501        CALL iom_put( 'iceshelf_cea', frcv(jpr_isf)%z3(:,:,1) ) 
     1502      ENDIF 
     1503 
     1504 
    14701505      IF( ln_mixcpl ) THEN 
    14711506         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:) 
     
    15041539      ENDIF 
    15051540 
     1541 
     1542      IF( srcv(jpr_icb)%laction )  THEN  
     1543         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1) 
     1544         rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runoffs 
     1545         CALL iom_put( 'iceberg_cea', frcv(jpr_icb)%z3(:,:,1) ) 
     1546      ENDIF 
     1547      IF( srcv(jpr_isf)%laction )  THEN 
     1548        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting)   
     1549        CALL iom_put( 'iceshelf_cea', frcv(jpr_isf)%z3(:,:,1) ) 
     1550      ENDIF 
     1551 
     1552 
    15061553      IF( ln_mixcpl ) THEN 
    15071554         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:) 
     
    15701617         IF( iom_use('hflx_cal_cea') )   CALL iom_put( 'hflx_cal_cea', - frcv(jpr_cal)%z3(:,:,1) * lfus )   ! heat flux from calving 
    15711618      ENDIF 
     1619 
     1620!!chris      
     1621!!    The heat content associated to the ice shelf in removed in the routine sbcisf.F90 
     1622      ! 
     1623      IF( srcv(jpr_icb)%laction )  zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting 
     1624      ! 
     1625!!      ! 
    15721626 
    15731627#if defined key_lim3       
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r7494 r7607  
    5555   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base 
    5656 
     57   LOGICAL, PUBLIC ::   l_isfcpl = .false.       ! isf recieved from oasis 
     58 
    5759 
    5860   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ? 
     
    9496    REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d 
    9597    REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d 
     98    REAL(wp)                            :: zhisf 
     99 
    96100      ! 
    97101      !!--------------------------------------------------------------------- 
     
    142146            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
    143147         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN 
    144             ALLOCATE( sf_rnfisf(1), STAT=ierror ) 
    145             ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 
    146             CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
     148            IF( .NOT.l_isfcpl ) THEN 
     149               ALLOCATE( sf_rnfisf(1), STAT=ierror ) 
     150               ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 
     151               CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
     152             ENDIF 
    147153 
    148154            !: read effective lenght (BG03) 
     
    185191             
    186192            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 
    187             ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 
    188             ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 
    189             ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 
    190             CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
    191             !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' ) 
     193            IF( .NOT.l_isfcpl ) THEN 
     194               ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 
     195               ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 
     196               ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 
     197               CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
     198               !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' ) 
     199            ENDIF 
    192200         END IF 
    193201          
     
    242250            CALL iom_put('vtbl',vtbl(:,:)) 
    243251            ! compute fwf and heat flux 
    244             CALL sbc_isf_cav (kt) 
     252            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt) 
     253            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * lfusisf              ! heat        flux 
     254            ENDIF 
    245255 
    246256         ELSE IF (nn_isf == 2) THEN 
     
    251261         ELSE IF (nn_isf == 3) THEN 
    252262            ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    253             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
    254             fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
     263            IF( .NOT.l_isfcpl ) THEN 
     264               CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
     265               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
     266            ENDIF 
    255267            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
    256268            stbl(:,:)   = soce 
     
    258270         ELSE IF (nn_isf == 4) THEN 
    259271            ! specified fwf and heat flux forcing beneath the ice shelf 
    260             CALL fld_read ( kt, nn_fsbc, sf_fwfisf   ) 
    261             !CALL fld_read ( kt, nn_fsbc, sf_qisf   ) 
    262             fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf 
     272            IF( .NOT.l_isfcpl ) THEN 
     273               CALL fld_read ( kt, nn_fsbc, sf_fwfisf   ) 
     274               !CALL fld_read ( kt, nn_fsbc, sf_qisf   ) 
     275               fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf 
     276            ENDIF 
    263277            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
    264278            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux 
     
    288302         CALL lbc_lnk(qisf(:,:)     ,'T',1.) 
    289303 
    290 !============================================================================================================================================= 
    291          IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 
     304         ! Diagnostics 
     305         IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 
     306            ! 
    292307            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 
    293308            CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        ) 
    294  
     309            ! 
    295310            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s) 
    296311            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2) 
    297312            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2) 
    298313            zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2) 
    299  
     314            ! 
    300315            DO jj = 1,jpj 
    301316               DO ji = 1,jpi 
     
    303318                  ikb = misfkb(ji,jj) 
    304319                  DO jk = ikt, ikb - 1 
    305                      zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
    306                      zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
    307                      zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
     320                     zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
     321                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf(ji,jj)    * zhisf 
     322                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf 
     323                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf(ji,jj)      * zhisf 
    308324                  END DO 
    309                   zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
    310                   zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
    311                   zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk) 
     325                  jk = ikb 
     326                  zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk) 
     327                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * zhisf * ralpha(ji,jj)  
     328                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf * ralpha(ji,jj) 
     329                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * zhisf * ralpha(ji,jj) 
    312330               END DO 
    313331            END DO 
    314  
    315             CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:)) 
    316             CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:)) 
    317             CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:)) 
    318             CALL iom_put('qhcisf'   , zqhcisf2d (:,:  )) 
    319  
     332            ! 
     333            CALL iom_put( 'fwfisf3d' , zfwfisf3d (:,:,:) ) 
     334            CALL iom_put( 'qlatisf3d', zqlatisf3d(:,:,:) ) 
     335            CALL iom_put( 'qhcisf3d' , zqhcisf3d (:,:,:) ) 
     336            CALL iom_put( 'qhcisf'   , zqhcisf2d (:,:  ) ) 
     337            ! 
    320338            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 
    321339            CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        ) 
     340            ! 
    322341         END IF 
    323 !============================================================================================================================================= 
    324  
     342         ! 
    325343         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
    326344            IF( ln_rstart .AND.    &                     ! Restart: read in restart file 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/stpctl.F90

    r6204 r7607  
    180180      ENDIF 
    181181 
    182 9200  FORMAT('it:', i8, ' iter:', i4, ' r: ',e16.10, ' b: ',e16.10 ) 
    183 9300  FORMAT(' it :', i8, ' ssh2: ', e16.10, ' Umax: ',e16.10,' Smin: ',e16.10) 
     1829200  FORMAT('it:', i8, ' iter:', i4, ' r: ',d23.16, ' b: ',d23.16 ) 
     1839300  FORMAT(' it :', i8, ' ssh2: ', d23.16, ' Umax: ',d23.16,' Smin: ',d23.16) 
    184184      ! 
    185185   END SUBROUTINE stp_ctl 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90

    r6943 r7607  
    1111   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90 
    1212   !!                  !  2011-02  (J. Simeon, J.Orr ) update O2 solubility constants 
     13   !!             3.6  !  2016-03  (O. Aumont) Change chemistry to MOCSY standards 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_pisces 
    1516   !!---------------------------------------------------------------------- 
    16    !!   'key_pisces'                                       PISCES bio-model 
     17   !!   'key_pisces*'                                      PISCES bio-model 
    1718   !!---------------------------------------------------------------------- 
    1819   !!   p4z_che      :  Sea water chemistry computed following OCMIP protocol 
     
    2223   USE sms_pisces    !  PISCES Source Minus Sink variables 
    2324   USE lib_mpp       !  MPP library 
     25   USE eosbn2, ONLY : nn_eos 
    2426 
    2527   IMPLICIT NONE 
    2628   PRIVATE 
    2729 
    28    PUBLIC   p4z_che         ! 
    29    PUBLIC   p4z_che_alloc   ! 
    30  
    31    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sio3eq   ! chemistry of Si 
    32    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 
     30   PUBLIC   p4z_che          ! 
     31   PUBLIC   p4z_che_alloc    ! 
     32   PUBLIC   p4z_che_ahini    ! 
     33   PUBLIC   p4z_che_solve_hi ! 
     34 
     35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: sio3eq   ! chemistry of Si 
     36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: fekeq    ! chemistry of Fe 
     37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: chemc    ! Solubilities of O2 and CO2 
     38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: chemo2    ! Solubilities of O2 and CO2 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: fesol    ! solubility of Fe 
    3540   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature 
     41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   salinprac  ! Practical salinity 
     42 
     43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ??? 
     44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akw3       !: ??? 
     45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akf3       !: ??? 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aks3       !: ??? 
     47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak1p3      !: ??? 
     48   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak2p3      !: ??? 
     49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak3p3      !: ??? 
     50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksi3      !: ??? 
     51   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   borat      !: ??? 
     52   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fluorid    !: ??? 
     53   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sulfat     !: ??? 
     54 
     55   !!* Variable for chemistry of the CO2 cycle 
    3656 
    3757   REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm 
    3858 
    39    REAL(wp) ::   salchl = 1. / 1.80655    ! conversion factor for salinity --> chlorinity (Wooster et al. 1969) 
    4059   REAL(wp) ::   o2atm  = 1. / ( 1000. * 0.20946 )   
    4160 
    42    REAL(wp) ::   rgas   = 83.14472       ! universal gas constants 
    43    REAL(wp) ::   oxyco  = 1. / 22.4144   ! converts from liters of an ideal gas to moles 
    44  
    45    REAL(wp) ::   bor1   = 0.00023        ! borat constants 
    46    REAL(wp) ::   bor2   = 1. / 10.82 
    47  
    48    REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate 
    49    REAL(wp) ::   st2    =  1./96.062    !  (Morris & Riley 1966) 
    50  
    51    REAL(wp) ::   ft1    =    0.000067   ! constants for calculate concentrations for fluorides 
    52    REAL(wp) ::   ft2    = 1./18.9984    ! (Dickson & Riley 1979 ) 
    53  
    54    !                                    ! volumetric solubility constants for o2 in ml/L   
    55    REAL(wp) ::   ox0    =  2.00856      ! from Table 1 for Eq 8 of Garcia and Gordon, 1992. 
    56    REAL(wp) ::   ox1    =  3.22400      ! corrects for moisture and fugacity, but not total atmospheric pressure 
    57    REAL(wp) ::   ox2    =  3.99063      !      Original PISCES code noted this was a solubility, but  
    58    REAL(wp) ::   ox3    =  4.80299      ! was in fact a bunsen coefficient with units L-O2/(Lsw atm-O2) 
    59    REAL(wp) ::   ox4    =  9.78188e-1   ! Hence, need to divide EXP( zoxy ) by 1000, ml-O2 => L-O2 
    60    REAL(wp) ::   ox5    =  1.71069      ! and atcox = 0.20946 to add the 1/atm dimension. 
    61    REAL(wp) ::   ox6    = -6.24097e-3    
    62    REAL(wp) ::   ox7    = -6.93498e-3  
    63    REAL(wp) ::   ox8    = -6.90358e-3 
    64    REAL(wp) ::   ox9    = -4.29155e-3  
    65    REAL(wp) ::   ox10   = -3.11680e-7  
     61   REAL(wp) ::   rgas   = 83.14472      ! universal gas constants 
     62   REAL(wp) ::   oxyco  = 1. / 22.4144  ! converts from liters of an ideal gas to moles 
    6663 
    6764   !                                    ! coeff. for seawater pressure correction : millero 95 
    6865   !                                    ! AGRIF doesn't like the DATA instruction 
    69    REAL(wp) :: devk11  = -25.5 
    70    REAL(wp) :: devk12  = -15.82 
    71    REAL(wp) :: devk13  = -29.48 
    72    REAL(wp) :: devk14  = -25.60 
    73    REAL(wp) :: devk15  = -48.76 
     66   REAL(wp) :: devk10  = -25.5 
     67   REAL(wp) :: devk11  = -15.82 
     68   REAL(wp) :: devk12  = -29.48 
     69   REAL(wp) :: devk13  = -20.02 
     70   REAL(wp) :: devk14  = -18.03 
     71   REAL(wp) :: devk15  = -9.78 
     72   REAL(wp) :: devk16  = -48.76 
     73   REAL(wp) :: devk17  = -14.51 
     74   REAL(wp) :: devk18  = -23.12 
     75   REAL(wp) :: devk19  = -26.57 
     76   REAL(wp) :: devk110  = -29.48 
    7477   ! 
    75    REAL(wp) :: devk21  = 0.1271 
    76    REAL(wp) :: devk22  = -0.0219 
    77    REAL(wp) :: devk23  = 0.1622 
    78    REAL(wp) :: devk24  = 0.2324 
    79    REAL(wp) :: devk25  = 0.5304 
     78   REAL(wp) :: devk20  = 0.1271 
     79   REAL(wp) :: devk21  = -0.0219 
     80   REAL(wp) :: devk22  = 0.1622 
     81   REAL(wp) :: devk23  = 0.1119 
     82   REAL(wp) :: devk24  = 0.0466 
     83   REAL(wp) :: devk25  = -0.0090 
     84   REAL(wp) :: devk26  = 0.5304 
     85   REAL(wp) :: devk27  = 0.1211 
     86   REAL(wp) :: devk28  = 0.1758 
     87   REAL(wp) :: devk29  = 0.2020 
     88   REAL(wp) :: devk210  = 0.1622 
    8089   ! 
     90   REAL(wp) :: devk30  = 0. 
    8191   REAL(wp) :: devk31  = 0. 
    82    REAL(wp) :: devk32  = 0. 
    83    REAL(wp) :: devk33  = 2.608E-3 
    84    REAL(wp) :: devk34  = -3.6246E-3 
    85    REAL(wp) :: devk35  = 0. 
     92   REAL(wp) :: devk32  = 2.608E-3 
     93   REAL(wp) :: devk33  = -1.409e-3 
     94   REAL(wp) :: devk34  = 0.316e-3 
     95   REAL(wp) :: devk35  = -0.942e-3 
     96   REAL(wp) :: devk36  = 0. 
     97   REAL(wp) :: devk37  = -0.321e-3 
     98   REAL(wp) :: devk38  = -2.647e-3 
     99   REAL(wp) :: devk39  = -3.042e-3 
     100   REAL(wp) :: devk310  = -2.6080e-3 
    86101   ! 
    87    REAL(wp) :: devk41  = -3.08E-3 
    88    REAL(wp) :: devk42  = 1.13E-3 
    89    REAL(wp) :: devk43  = -2.84E-3 
    90    REAL(wp) :: devk44  = -5.13E-3 
    91    REAL(wp) :: devk45  = -11.76E-3 
     102   REAL(wp) :: devk40  = -3.08E-3 
     103   REAL(wp) :: devk41  = 1.13E-3 
     104   REAL(wp) :: devk42  = -2.84E-3 
     105   REAL(wp) :: devk43  = -5.13E-3 
     106   REAL(wp) :: devk44  = -4.53e-3 
     107   REAL(wp) :: devk45  = -3.91e-3 
     108   REAL(wp) :: devk46  = -11.76e-3 
     109   REAL(wp) :: devk47  = -2.67e-3 
     110   REAL(wp) :: devk48  = -5.15e-3 
     111   REAL(wp) :: devk49  = -4.08e-3 
     112   REAL(wp) :: devk410  = -2.84e-3 
    92113   ! 
    93    REAL(wp) :: devk51  = 0.0877E-3 
    94    REAL(wp) :: devk52  = -0.1475E-3      
    95    REAL(wp) :: devk53  = 0. 
    96    REAL(wp) :: devk54  = 0.0794E-3       
    97    REAL(wp) :: devk55  = 0.3692E-3       
     114   REAL(wp) :: devk50  = 0.0877E-3 
     115   REAL(wp) :: devk51  = -0.1475E-3      
     116   REAL(wp) :: devk52  = 0. 
     117   REAL(wp) :: devk53  = 0.0794E-3       
     118   REAL(wp) :: devk54  = 0.09e-3 
     119   REAL(wp) :: devk55  = 0.054e-3 
     120   REAL(wp) :: devk56  = 0.3692E-3 
     121   REAL(wp) :: devk57  = 0.0427e-3 
     122   REAL(wp) :: devk58  = 0.09e-3 
     123   REAL(wp) :: devk59  = 0.0714e-3 
     124   REAL(wp) :: devk510  = 0.0 
     125   ! 
     126   ! General parameters 
     127   REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp 
     128   REAL(wp), PARAMETER :: pp_ln10 = 2.302585092994045684018_wp 
     129 
     130   ! Maximum number of iterations for each method 
     131   INTEGER, PARAMETER :: jp_maxniter_atgen    = 20 
     132 
     133   ! Bookkeeping variables for each method 
     134   ! - SOLVE_AT_GENERAL 
     135   INTEGER :: niter_atgen    = jp_maxniter_atgen 
    98136 
    99137   !!* Substitution 
     
    115153      !!--------------------------------------------------------------------- 
    116154      INTEGER  ::   ji, jj, jk 
    117       REAL(wp) ::   ztkel, zt   , zt2  , zsal  , zsal2 , zbuf1 , zbuf2 
     155      REAL(wp) ::   ztkel, ztkel1, zt , zsal  , zsal2 , zbuf1 , zbuf2 
    118156      REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5 
    119157      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2 
    120158      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat 
    121       REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1  , za2 
     159      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1, za2 
    122160      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw 
     161      REAL(wp) ::   zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi 
    123162      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1 
     163      REAL(wp) ::   total2free, free2SWS, total2SWS, SWS2total 
     164 
    124165      !!--------------------------------------------------------------------- 
    125166      ! 
    126167      IF( nn_timing == 1 )  CALL timing_start('p4z_che') 
     168      ! 
     169      ! Computation of chemical constants require practical salinity 
     170      ! Thus, when TEOS08 is used, absolute salinity is converted to  
     171      ! practical salinity 
     172      ! ------------------------------------------------------------- 
     173      IF (nn_eos == -1) THEN 
     174         salinprac(:,:,:) = tsn(:,:,:,jp_sal) * 35.0 / 35.16504 
     175      ELSE 
     176         salinprac(:,:,:) = tsn(:,:,:,jp_sal) 
     177      ENDIF 
     178 
    127179      ! 
    128180      ! Computations of chemical constants require in situ temperature 
     
    135187            DO ji = 1, jpi 
    136188               zpres = fsdept(ji,jj,jk) / 1000. 
    137                za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (tsn(ji,jj,jk,jp_sal) - 35.0) ) 
     189               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (salinprac(ji,jj,jk) - 35.0) ) 
    138190               za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 ) 
    139191               tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2 
     
    151203            ztkel = tempis(ji,jj,1) + 273.15 
    152204            zt    = ztkel * 0.01 
    153             zt2   = zt * zt 
    154             zsal  = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35. 
    155             zsal2 = zsal * zsal 
    156             zlogt = LOG( zt ) 
     205            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35. 
    157206            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 
    158207            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
    159208            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    & 
    160209            &       + 0.0047036e-4*ztkel**2) 
    161             !                             ! SET SOLUBILITIES OF O2 AND CO2  
    162             chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000. ! mol/(kg uatm) 
     210            chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm) 
    163211            chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 
    164212            chemc(ji,jj,3) = 57.7 - 0.118*ztkel 
     
    176224            DO ji = 1, jpi 
    177225              ztkel = tempis(ji,jj,jk) + 273.15 
    178               zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35. 
     226              zsal  = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35. 
    179227              zsal2 = zsal * zsal 
    180228              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
     
    183231              ztgg4 = ztgg3 * ztgg 
    184232              ztgg5 = ztgg4 * ztgg 
    185               zoxy  = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5   & 
    186                      + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) +  ox10 * zsal2 
     233 
     234              zoxy  = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3    & 
     235              &       + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3   & 
     236              &       - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 )   & 
     237              &       - 3.11680e-7 * zsal2 
    187238              chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox     ! mol/(L atm) 
    188239            END DO 
     
    209260               ! SET ABSOLUTE TEMPERATURE 
    210261               ztkel   = tempis(ji,jj,jk) + 273.15 
    211                zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35. 
     262               zsal    = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35. 
    212263               zsqrt  = SQRT( zsal ) 
    213264               zsal15  = zsqrt * zsal 
     
    220271 
    221272               ! CHLORINITY (WOOSTER ET AL., 1969) 
    222                zcl     = zsal * salchl 
     273               zcl     = zsal / 1.80655 
    223274 
    224275               ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 
    225                zst     = st1 * zcl * st2 
     276               zst     = 0.14 * zcl /96.062 
    226277 
    227278               ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 
    228                zft     = ft1 * zcl * ft2 
     279               zft     = 0.000067 * zcl /18.9984 
    229280 
    230281               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 
     
    234285               &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         & 
    235286               &         + LOG(1.0 - 0.001005 * zsal)) 
    236                ! 
    237                aphscale(ji,jj,jk) = ( 1. + zst / zcks ) 
    238287 
    239288               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 
     
    249298               &      * zlogt + 0.053105*zsqrt*ztkel 
    250299 
    251  
    252300               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO  
    253301               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 
     
    257305                  - 0.01781*zsal + 0.0001122*zsal*zsal) 
    258306 
    259                ! PKW (H2O) (DICKSON AND RILEY, 1979) 
    260                zckw = -13847.26*ztr + 148.9652 - 23.6521 * zlogt    &  
    261                &     + (118.67*ztr - 5.977 + 1.0495 * zlogt)        & 
    262                &     * zsqrt - 0.01615 * zsal 
     307               ! PKW (H2O) (MILLERO, 1995) from composite data 
     308               zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    & 
     309                         - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal 
     310 
     311               ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995) 
     312              zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   & 
     313              &          + (-106.736*ztr + 0.69171) * zsqrt       & 
     314              &          + (-0.65643*ztr - 0.01844) * zsal 
     315 
     316              zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  & 
     317              &          + (-160.340*ztr + 1.3566)*zsqrt          & 
     318              &          + (0.37335*ztr - 0.05778)*zsal 
     319 
     320              zck3p    = -3070.75*ztr - 18.126                    & 
     321              &          + (17.27039*ztr + 2.81197) * zsqrt       & 
     322              &          + (-44.99486*ztr - 0.09984) * zsal  
     323 
     324              ! CONSTANT FOR SILICATE, MILLERO (1995) 
     325              zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   & 
     326              &          + (-458.79*ztr + 3.5913) * zisqrt       & 
     327              &          + (188.74*ztr - 1.5998) * zis           & 
     328              &          + (-12.1652*ztr + 0.07871) * zis2       & 
     329              &          + LOG(1.0 - 0.001005*zsal) 
    263330 
    264331               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
     
    268335                  &      - 0.07711*zsal + 0.0041249*zsal15 
    269336 
     337               ! CONVERT FROM DIFFERENT PH SCALES 
     338               total2free  = 1.0/(1.0 + zst/zcks) 
     339               free2SWS    = 1. + zst/zcks + zft/(zckf*total2free) 
     340               total2SWS   = total2free * free2SWS 
     341               SWS2total   = 1.0 / total2SWS 
     342 
    270343               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 
    271                zak1    = 10**(zck1) 
    272                zak2    = 10**(zck2) 
    273                zakb    = EXP( zckb  ) 
     344               zak1    = 10**(zck1) * total2SWS 
     345               zak2    = 10**(zck2) * total2SWS 
     346               zakb    = EXP( zckb ) * total2SWS 
    274347               zakw    = EXP( zckw ) 
    275348               zaksp1  = 10**(zaksp0) 
     349               zak1p   = exp( zck1p ) 
     350               zak2p   = exp( zck2p ) 
     351               zak3p   = exp( zck3p ) 
     352               zaksi   = exp( zcksi ) 
     353               zckf    = zckf * total2SWS 
    276354 
    277355               ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) 
     
    285363               !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 
    286364               !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 
    287                zcpexp  = zpres /(rgas*ztkel) 
    288                zcpexp2 = zpres * zpres/(rgas*ztkel) 
     365               zcpexp  = zpres / (rgas*ztkel) 
     366               zcpexp2 = zpres * zcpexp 
    289367 
    290368               ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE 
     
    292370               !        (CF. BROECKER ET AL., 1982) 
    293371 
    294                zbuf1  = -     ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 
     372               zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc ) 
     373               zbuf2  = 0.5 * ( devk40 + devk50 * ztc ) 
     374               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     375 
     376               zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 
    295377               zbuf2  = 0.5 * ( devk41 + devk51 * ztc ) 
    296                ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     378               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
    297379 
    298380               zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc ) 
    299381               zbuf2  = 0.5 * ( devk42 + devk52 * ztc ) 
    300                ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     382               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
    301383 
    302384               zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc ) 
    303385               zbuf2  = 0.5 * ( devk43 + devk53 * ztc ) 
    304                akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     386               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
    305387 
    306388               zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc ) 
    307389               zbuf2  = 0.5 * ( devk44 + devk54 * ztc ) 
    308                akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
    309  
     390               aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     391 
     392               zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 
     393               zbuf2  = 0.5 * ( devk45 + devk55 * ztc ) 
     394               akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     395 
     396               zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc ) 
     397               zbuf2  = 0.5 * ( devk47 + devk57 * ztc ) 
     398               ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     399 
     400               zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc ) 
     401               zbuf2  = 0.5 * ( devk48 + devk58 * ztc ) 
     402               ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     403 
     404               zbuf1  =     - ( devk19 + devk29 * ztc + devk39 * ztc * ztc ) 
     405               zbuf2  = 0.5 * ( devk49 + devk59 * ztc ) 
     406               ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     407 
     408               zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc ) 
     409               zbuf2  = 0.5 * ( devk410 + devk510 * ztc ) 
     410               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
     411 
     412               ! CONVERT FROM DIFFERENT PH SCALES 
     413               total2free  = 1.0/(1.0 + zst/aks3(ji,jj,jk)) 
     414               free2SWS    = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk) 
     415               total2SWS   = total2free * free2SWS 
     416               SWS2total   = 1.0 / total2SWS 
     417 
     418               ! Convert to total scale 
     419               ak13(ji,jj,jk)  = ak13(ji,jj,jk)  * SWS2total 
     420               ak23(ji,jj,jk)  = ak23(ji,jj,jk)  * SWS2total 
     421               akb3(ji,jj,jk)  = akb3(ji,jj,jk)  * SWS2total 
     422               akw3(ji,jj,jk)  = akw3(ji,jj,jk)  * SWS2total 
     423               ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total 
     424               ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total 
     425               ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total 
     426               aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total 
     427               akf3(ji,jj,jk)  = akf3(ji,jj,jk)  / total2free 
    310428 
    311429               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE  
    312430               !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO 
    313431               !        (P. 1285) AND BERNER (1976) 
    314                zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 
    315                zbuf2  = 0.5 * ( devk45 + devk55 * ztc ) 
     432               zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc ) 
     433               zbuf2  = 0.5 * ( devk46 + devk56 * ztc ) 
    316434               aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 
    317435 
    318                ! TOTAL BORATE CONCENTR. [MOLES/L] 
    319                borat(ji,jj,jk) = bor1 * zcl * bor2 
     436               ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L] 
     437               borat(ji,jj,jk) = 0.0002414 * zcl / 10.811 
     438               sulfat(ji,jj,jk) = zst 
     439               fluorid(ji,jj,jk) = zft  
    320440 
    321441               ! Iron and SIO3 saturation concentration from ... 
    322442               sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6 
    323                fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ( 273.15 + ztc ) ) 
    324  
     443               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel ) 
     444 
     445               ! Liu and Millero (1999) only valid 5 - 50 degC 
     446               ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16 
     447               fesol(ji,jj,jk,1) = 10**((-13.486) - (0.1856* (zis**0.5)) + (0.3073*zis) + (5254.0/ztkel1)) 
     448               fesol(ji,jj,jk,2) = 10**(2.517 - (0.885*(zis**0.5)) + (0.2139 * zis) - (1320.0/ztkel1) ) 
     449               fesol(ji,jj,jk,3) = 10**(0.4511 - (0.3305*(zis**0.5)) - (1996.0/ztkel1) ) 
     450               fesol(ji,jj,jk,4) = 10**(-0.2965 - (0.7881*(zis**0.5)) - (4086.0/ztkel1) ) 
     451               fesol(ji,jj,jk,5) = 10**(4.4466 - (0.8505*(zis**0.5)) - (7980.0/ztkel1) ) 
    325452            END DO 
    326453         END DO 
     
    331458   END SUBROUTINE p4z_che 
    332459 
     460   SUBROUTINE p4z_che_ahini( p_hini ) 
     461      !!--------------------------------------------------------------------- 
     462      !!                     ***  ROUTINE ahini_for_at  *** 
     463      !! 
     464      !! Subroutine returns the root for the 2nd order approximation of the 
     465      !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic  
     466      !! polynomial) around the local minimum, if it exists. 
     467      !! Returns * 1E-03_wp if p_alkcb <= 0 
     468      !!         * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot 
     469      !!         * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot 
     470      !!                    and the 2nd order approximation does not have  
     471      !!                    a solution 
     472      !!--------------------------------------------------------------------- 
     473      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  ::  p_hini 
     474      INTEGER  ::   ji, jj, jk 
     475      REAL(wp)  ::  zca1, zba1 
     476      REAL(wp)  ::  zd, zsqrtd, zhmin 
     477      REAL(wp)  ::  za2, za1, za0 
     478      REAL(wp)  ::  p_dictot, p_bortot, p_alkcb  
     479 
     480      IF( nn_timing == 1 )  CALL timing_start('p4z_che_ahini') 
     481      ! 
     482      DO jk = 1, jpk 
     483        DO jj = 1, jpj 
     484          DO ji = 1, jpi 
     485            p_alkcb  = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 
     486            p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn) 
     487            p_bortot = borat(ji,jj,jk) 
     488            IF (p_alkcb <= 0.) THEN 
     489                p_hini(ji,jj,jk) = 1.e-3 
     490            ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN 
     491                p_hini(ji,jj,jk) = 1.e-10_wp 
     492            ELSE 
     493                zca1 = p_dictot/( p_alkcb + rtrn ) 
     494                zba1 = p_bortot/ (p_alkcb + rtrn ) 
     495           ! Coefficients of the cubic polynomial 
     496                za2 = aKb3(ji,jj,jk)*(1. - zba1) + ak13(ji,jj,jk)*(1.-zca1) 
     497                za1 = ak13(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - zca1)    & 
     498                &     + ak13(ji,jj,jk)*ak23(ji,jj,jk)*(1. - (zca1+zca1)) 
     499                za0 = ak13(ji,jj,jk)*ak23(ji,jj,jk)*akb3(ji,jj,jk)*(1. - zba1 - (zca1+zca1)) 
     500                                        ! Taylor expansion around the minimum 
     501                zd = za2*za2 - 3.*za1   ! Discriminant of the quadratic equation 
     502                                        ! for the minimum close to the root 
     503 
     504                IF(zd > 0.) THEN        ! If the discriminant is positive 
     505                  zsqrtd = SQRT(zd) 
     506                  IF(za2 < 0) THEN 
     507                    zhmin = (-za2 + zsqrtd)/3. 
     508                  ELSE 
     509                    zhmin = -za1/(za2 + zsqrtd) 
     510                  ENDIF 
     511                  p_hini(ji,jj,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd) 
     512                ELSE 
     513                  p_hini(ji,jj,jk) = 1.e-7 
     514                ENDIF 
     515             ! 
     516             ENDIF 
     517          END DO 
     518        END DO 
     519      END DO 
     520      ! 
     521      IF( nn_timing == 1 )  CALL timing_stop('p4z_che_ahini') 
     522      ! 
     523   END SUBROUTINE p4z_che_ahini 
     524 
     525   !=============================================================================== 
     526   SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup ) 
     527 
     528   ! Subroutine returns the lower and upper bounds of "non-water-selfionization" 
     529   ! contributions to total alkalinity (the infimum and the supremum), i.e 
     530   ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+]) 
     531 
     532   ! Argument variables 
     533   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf 
     534   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup 
     535 
     536   p_alknw_inf(:,:,:) =  -trb(:,:,:,jppo4) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:)  & 
     537   &              - fluorid(:,:,:) 
     538   p_alknw_sup(:,:,:) =   (2. * trb(:,:,:,jpdic) + 2. * trb(:,:,:,jppo4) + trb(:,:,:,jpsil) )    & 
     539   &               * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:)  
     540 
     541   END SUBROUTINE anw_infsup 
     542 
     543 
     544   SUBROUTINE p4z_che_solve_hi( p_hini, zhi ) 
     545 
     546   ! Universal pH solver that converges from any given initial value, 
     547   ! determines upper an lower bounds for the solution if required 
     548 
     549   ! Argument variables 
     550   !-------------------- 
     551   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN)   :: p_hini 
     552   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  :: zhi 
     553 
     554   ! Local variables 
     555   !----------------- 
     556   INTEGER   ::  ji, jj, jk, jn 
     557   REAL(wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor 
     558   REAL(wp)  ::  zdelta, zh_delta 
     559   REAL(wp)  ::  zeqn, zdeqndh, zalka 
     560   REAL(wp)  ::  aphscale 
     561   REAL(wp)  ::  znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic 
     562   REAL(wp)  ::  znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor 
     563   REAL(wp)  ::  znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4 
     564   REAL(wp)  ::  znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil 
     565   REAL(wp)  ::  znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 
     566   REAL(wp)  ::  znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu 
     567   REAL(wp)  ::  zalk_wat, zdalk_wat 
     568   REAL(wp)  ::  zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit 
     569   LOGICAL   ::  l_exitnow 
     570   REAL(wp), PARAMETER :: pz_exp_threshold = 1.0 
     571   REAL(wp), POINTER, DIMENSION(:,:,:) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin 
     572 
     573   IF( nn_timing == 1 )  CALL timing_start('p4z_che_solve_hi') 
     574      !  Allocate temporary workspace 
     575   CALL wrk_alloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask ) 
     576   CALL wrk_alloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin ) 
     577 
     578   CALL anw_infsup( zalknw_inf, zalknw_sup ) 
     579 
     580   rmask(:,:,:) = tmask(:,:,:) 
     581   zhi(:,:,:)   = 0. 
     582 
     583   ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree 
     584   DO jk = 1, jpk 
     585      DO jj = 1, jpj 
     586         DO ji = 1, jpi 
     587            IF (rmask(ji,jj,jk) == 1.) THEN 
     588               p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn) 
     589               aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 
     590               zh_ini = p_hini(ji,jj,jk) 
     591 
     592               zdelta = (p_alktot-zalknw_inf(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 
     593 
     594               IF(p_alktot >= zalknw_inf(ji,jj,jk)) THEN 
     595                 zh_min(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_inf(ji,jj,jk) + SQRT(zdelta) ) 
     596               ELSE 
     597                 zh_min(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_inf(ji,jj,jk)) + SQRT(zdelta) ) / 2. 
     598               ENDIF 
     599 
     600               zdelta = (p_alktot-zalknw_sup(ji,jj,jk))**2 + 4.*akw3(ji,jj,jk)/aphscale 
     601 
     602               IF(p_alktot <= zalknw_sup(ji,jj,jk)) THEN 
     603                 zh_max(ji,jj,jk) = aphscale*(-(p_alktot-zalknw_sup(ji,jj,jk)) + SQRT(zdelta) ) / 2. 
     604               ELSE 
     605                 zh_max(ji,jj,jk) = 2.*akw3(ji,jj,jk) /( p_alktot-zalknw_sup(ji,jj,jk) + SQRT(zdelta) ) 
     606               ENDIF 
     607 
     608               zhi(ji,jj,jk) = MAX(MIN(zh_max(ji,jj,jk), zh_ini), zh_min(ji,jj,jk)) 
     609            ENDIF 
     610         END DO 
     611      END DO 
     612   END DO 
     613 
     614   zeqn_absmin(:,:,:) = HUGE(1._wp) 
     615 
     616   DO jn = 1, jp_maxniter_atgen  
     617   DO jk = 1, jpk 
     618      DO jj = 1, jpj 
     619         DO ji = 1, jpi 
     620            IF (rmask(ji,jj,jk) == 1.) THEN 
     621               zfact = rhop(ji,jj,jk) / 1000. + rtrn 
     622               p_alktot = trb(ji,jj,jk,jptal) / zfact 
     623               zdic  = trb(ji,jj,jk,jpdic) / zfact 
     624               zbot  = borat(ji,jj,jk) 
     625               zpt = trb(ji,jj,jk,jppo4) / zfact * po4r 
     626               zsit = trb(ji,jj,jk,jpsil) / zfact 
     627               zst = sulfat (ji,jj,jk) 
     628               zft = fluorid(ji,jj,jk) 
     629               aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk) 
     630               zh = zhi(ji,jj,jk) 
     631               zh_prev = zh 
     632 
     633               ! H2CO3 - HCO3 - CO3 : n=2, m=0 
     634               znumer_dic = 2.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk) 
     635               zdenom_dic = ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*(ak13(ji,jj,jk) + zh) 
     636               zalk_dic   = zdic * (znumer_dic/zdenom_dic) 
     637               zdnumer_dic = ak13(ji,jj,jk)*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh     & 
     638                             *(4.*ak13(ji,jj,jk)*ak23(ji,jj,jk) + zh*ak13(ji,jj,jk)) 
     639               zdalk_dic   = -zdic*(zdnumer_dic/zdenom_dic**2) 
     640 
     641 
     642               ! B(OH)3 - B(OH)4 : n=1, m=0 
     643               znumer_bor = akb3(ji,jj,jk) 
     644               zdenom_bor = akb3(ji,jj,jk) + zh 
     645               zalk_bor   = zbot * (znumer_bor/zdenom_bor) 
     646               zdnumer_bor = akb3(ji,jj,jk) 
     647               zdalk_bor   = -zbot*(zdnumer_bor/zdenom_bor**2) 
     648 
     649 
     650               ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1 
     651               znumer_po4 = 3.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)  & 
     652               &            + zh*(2.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh* ak1p3(ji,jj,jk)) 
     653               zdenom_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)     & 
     654               &            + zh*( ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh*(ak1p3(ji,jj,jk) + zh)) 
     655               zalk_po4   = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1 
     656               zdnumer_po4 = ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)  & 
     657               &             + zh*(4.*ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)         & 
     658               &             + zh*(9.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)*ak3p3(ji,jj,jk)                         & 
     659               &             + ak1p3(ji,jj,jk)*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk)                                & 
     660               &             + zh*(4.*ak1p3(ji,jj,jk)*ak2p3(ji,jj,jk) + zh * ak1p3(ji,jj,jk) ) ) ) 
     661               zdalk_po4   = -zpt * (zdnumer_po4/zdenom_po4**2) 
     662 
     663               ! H4SiO4 - H3SiO4 : n=1, m=0 
     664               znumer_sil = aksi3(ji,jj,jk) 
     665               zdenom_sil = aksi3(ji,jj,jk) + zh 
     666               zalk_sil   = zsit * (znumer_sil/zdenom_sil) 
     667               zdnumer_sil = aksi3(ji,jj,jk) 
     668               zdalk_sil   = -zsit * (zdnumer_sil/zdenom_sil**2) 
     669 
     670               ! HSO4 - SO4 : n=1, m=1 
     671               aphscale = 1.0 + zst/aks3(ji,jj,jk) 
     672               znumer_so4 = aks3(ji,jj,jk) * aphscale 
     673               zdenom_so4 = aks3(ji,jj,jk) * aphscale + zh 
     674               zalk_so4   = zst * (znumer_so4/zdenom_so4 - 1.) 
     675               zdnumer_so4 = aks3(ji,jj,jk) 
     676               zdalk_so4   = -zst * (zdnumer_so4/zdenom_so4**2) 
     677 
     678               ! HF - F : n=1, m=1 
     679               znumer_flu =  akf3(ji,jj,jk) 
     680               zdenom_flu =  akf3(ji,jj,jk) + zh 
     681               zalk_flu   =  zft * (znumer_flu/zdenom_flu - 1.) 
     682               zdnumer_flu = akf3(ji,jj,jk) 
     683               zdalk_flu   = -zft * (zdnumer_flu/zdenom_flu**2) 
     684 
     685               ! H2O - OH 
     686               aphscale = 1.0 + zst/aks3(ji,jj,jk) 
     687               zalk_wat   = akw3(ji,jj,jk)/zh - zh/aphscale 
     688               zdalk_wat  = -akw3(ji,jj,jk)/zh**2 - 1./aphscale 
     689 
     690               ! CALCULATE [ALK]([CO3--], [HCO3-]) 
     691               zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil   & 
     692               &      + zalk_so4 + zalk_flu                       & 
     693               &      + zalk_wat - p_alktot 
     694 
     695               zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil   & 
     696               &       + zalk_so4 + zalk_flu + zalk_wat) 
     697 
     698               zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 
     699               &         + zdalk_so4 + zdalk_flu + zdalk_wat 
     700 
     701               ! Adapt bracketing interval 
     702               IF(zeqn > 0._wp) THEN 
     703                 zh_min(ji,jj,jk) = zh_prev 
     704               ELSEIF(zeqn < 0._wp) THEN 
     705                 zh_max(ji,jj,jk) = zh_prev 
     706               ENDIF 
     707 
     708               IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jj,jk)) THEN 
     709               ! if the function evaluation at the current point is 
     710               ! not decreasing faster than with a bisection step (at least linearly) 
     711               ! in absolute value take one bisection step on [ph_min, ph_max] 
     712               ! ph_new = (ph_min + ph_max)/2d0 
     713               ! 
     714               ! In terms of [H]_new: 
     715               ! [H]_new = 10**(-ph_new) 
     716               !         = 10**(-(ph_min + ph_max)/2d0) 
     717               !         = SQRT(10**(-(ph_min + phmax))) 
     718               !         = SQRT(zh_max * zh_min) 
     719                  zh = SQRT(zh_max(ji,jj,jk) * zh_min(ji,jj,jk)) 
     720                  zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below 
     721               ELSE 
     722               ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH 
     723               !           = -zdeqndh * LOG(10) * [H] 
     724               ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10)) 
     725               ! 
     726               ! pH_new = pH_old + \deltapH 
     727               ! 
     728               ! [H]_new = 10**(-pH_new) 
     729               !         = 10**(-pH_old - \Delta pH) 
     730               !         = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10))) 
     731               !         = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10))) 
     732               !         = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old)) 
     733 
     734                  zh_lnfactor = -zeqn/(zdeqndh*zh_prev) 
     735 
     736                  IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN 
     737                     zh          = zh_prev*EXP(zh_lnfactor) 
     738                  ELSE 
     739                     zh_delta    = zh_lnfactor*zh_prev 
     740                     zh          = zh_prev + zh_delta 
     741                  ENDIF 
     742 
     743                  IF( zh < zh_min(ji,jj,jk) ) THEN 
     744                     ! if [H]_new < [H]_min 
     745                     ! i.e., if ph_new > ph_max then 
     746                     ! take one bisection step on [ph_prev, ph_max] 
     747                     ! ph_new = (ph_prev + ph_max)/2d0 
     748                     ! In terms of [H]_new: 
     749                     ! [H]_new = 10**(-ph_new) 
     750                     !         = 10**(-(ph_prev + ph_max)/2d0) 
     751                     !         = SQRT(10**(-(ph_prev + phmax))) 
     752                     !         = SQRT([H]_old*10**(-ph_max)) 
     753                     !         = SQRT([H]_old * zh_min) 
     754                     zh                = SQRT(zh_prev * zh_min(ji,jj,jk)) 
     755                     zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below 
     756                  ENDIF 
     757 
     758                  IF( zh > zh_max(ji,jj,jk) ) THEN 
     759                     ! if [H]_new > [H]_max 
     760                     ! i.e., if ph_new < ph_min, then 
     761                     ! take one bisection step on [ph_min, ph_prev] 
     762                     ! ph_new = (ph_prev + ph_min)/2d0 
     763                     ! In terms of [H]_new: 
     764                     ! [H]_new = 10**(-ph_new) 
     765                     !         = 10**(-(ph_prev + ph_min)/2d0) 
     766                     !         = SQRT(10**(-(ph_prev + ph_min))) 
     767                     !         = SQRT([H]_old*10**(-ph_min)) 
     768                     !         = SQRT([H]_old * zhmax) 
     769                     zh                = SQRT(zh_prev * zh_max(ji,jj,jk)) 
     770                     zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below 
     771                  ENDIF 
     772               ENDIF 
     773 
     774               zeqn_absmin(ji,jj,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jj,jk)) 
     775 
     776               ! Stop iterations once |\delta{[H]}/[H]| < rdel 
     777               ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel 
     778               ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)| 
     779 
     780               ! Alternatively: 
     781               ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))| 
     782               !             ~ 1/LOG(10) * |\Delta [H]|/[H] 
     783               !             < 1/LOG(10) * rdel 
     784 
     785               ! Hence |zeqn/(zdeqndh*zh)| < rdel 
     786 
     787               ! rdel <-- pp_rdel_ah_target 
     788               l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target) 
     789 
     790               IF(l_exitnow) THEN  
     791                  rmask(ji,jj,jk) = 0. 
     792               ENDIF 
     793 
     794               zhi(ji,jj,jk) =  zh 
     795 
     796               IF(jn >= jp_maxniter_atgen) THEN 
     797                  zhi(ji,jj,jk) = -1._wp 
     798               ENDIF 
     799 
     800            ENDIF 
     801         END DO 
     802      END DO 
     803   END DO 
     804   END DO 
     805   ! 
     806   CALL wrk_dealloc( jpi, jpj, jpk, zalknw_inf, zalknw_sup, rmask ) 
     807   CALL wrk_dealloc( jpi, jpj, jpk, zh_min, zh_max, zeqn_absmin ) 
     808 
     809 
     810   IF( nn_timing == 1 )  CALL timing_stop('p4z_che_solve_hi') 
     811 
     812 
     813   END SUBROUTINE p4z_che_solve_hi 
    333814 
    334815   INTEGER FUNCTION p4z_che_alloc() 
     
    336817      !!                     ***  ROUTINE p4z_che_alloc  *** 
    337818      !!---------------------------------------------------------------------- 
    338       ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk),   & 
    339       &         tempis(jpi,jpj,jpk), STAT=p4z_che_alloc ) 
     819      INTEGER ::   ierr(3)        ! Local variables 
     820      !!---------------------------------------------------------------------- 
     821 
     822      ierr(:) = 0 
     823 
     824      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) ) 
     825 
     826      ALLOCATE( akb3(jpi,jpj,jpk)     , tempis(jpi, jpj, jpk),       & 
     827         &      akw3(jpi,jpj,jpk)     , borat (jpi,jpj,jpk)  ,       & 
     828         &      aks3(jpi,jpj,jpk)     , akf3(jpi,jpj,jpk)    ,       & 
     829         &      ak1p3(jpi,jpj,jpk)    , ak2p3(jpi,jpj,jpk)   ,       & 
     830         &      ak3p3(jpi,jpj,jpk)    , aksi3(jpi,jpj,jpk)   ,       & 
     831         &      fluorid(jpi,jpj,jpk)  , sulfat(jpi,jpj,jpk)  ,       & 
     832         &      salinprac(jpi,jpj,jpk),                 STAT=ierr(2) ) 
     833 
     834      ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) ) 
     835 
     836      !* Variable for chemistry of the CO2 cycle 
     837      p4z_che_alloc = MAXVAL( ierr ) 
    340838      ! 
    341839      IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') 
     
    355853 
    356854   !!====================================================================== 
    357 END MODULE p4zche 
     855END MODULE  p4zche 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90

    r6943 r7607  
    8787      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact 
    8888      REAL(wp) ::   zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff 
    89       REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 
     89      REAL(wp) ::   zph, zdic, zsch_o2, zsch_co2 
    9090      REAL(wp) ::   zyr_dec, zdco2dt 
    9191      CHARACTER (len=25) :: charout 
     
    122122#endif 
    123123 
    124       DO jm = 1, 10 
    125 !CDIR NOVERRCHK 
    126          DO jj = 1, jpj 
    127 !CDIR NOVERRCHK 
    128             DO ji = 1, jpi 
    129  
    130                ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 
    131                zbot  = borat(ji,jj,1) 
    132                zfact = rhop(ji,jj,1) / 1000. + rtrn 
    133                zdic  = trb(ji,jj,1,jpdic) / zfact 
    134                zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 
    135                zalka = trb(ji,jj,1,jptal) / zfact 
    136  
    137                ! CALCULATE [ALK]([CO3--], [HCO3-]) 
    138                zalk  = zalka - (  akw3(ji,jj,1) / zph - zph / aphscale(ji,jj,1)    & 
    139                &       + zbot / ( 1.+ zph / akb3(ji,jj,1) )  ) 
    140  
    141                ! CALCULATE [H+] AND [H2CO3] 
    142                zah2   = SQRT(  (zdic-zalk)**2 + 4.* ( zalk * ak23(ji,jj,1)   & 
    143                   &                                        / ak13(ji,jj,1) ) * ( 2.* zdic - zalk )  ) 
    144                zah2   = 0.5 * ak13(ji,jj,1) / zalk * ( ( zdic - zalk ) + zah2 ) 
    145                zh2co3(ji,jj) = ( 2.* zdic - zalk ) / ( 2.+ ak13(ji,jj,1) / zah2 ) * zfact 
    146                hi(ji,jj,1)   = zah2 * zfact 
    147             END DO 
     124      DO jj = 1, jpj 
     125         DO ji = 1, jpi 
     126            ! DUMMY VARIABLES FOR DIC, H+, AND BORATE 
     127            zfact = rhop(ji,jj,1) / 1000. + rtrn 
     128            zdic  = trb(ji,jj,1,jpdic) 
     129            zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact 
     130            ! CALCULATE [H2CO3] 
     131            zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2) 
    148132         END DO 
    149133      END DO 
    150  
    151134 
    152135      ! -------------- 
     
    254237      ENDIF 
    255238      ! 
     239#if defined key_cpl_carbon_cycle 
     240      ! change units for carbon cycle coupling 
     241      oce_co2(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r ! in molC/m2/s 
     242#endif 
     243      ! 
    256244      CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx, zpco2atm ) 
    257245      ! 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

    r6943 r7607  
    2323   USE sms_pisces      !  PISCES Source Minus Sink variables 
    2424   USE prtctl_trc      !  print control for debugging 
     25   USE p4zche          !  Chemical model 
    2526   USE iom             !  I/O manager 
    2627 
     
    6061      ! 
    6162      INTEGER, INTENT(in) ::   kt, knt ! ocean time step 
    62       INTEGER  ::   ji, jj, jk, jn 
    63       REAL(wp) ::   zalk, zdic, zph, zah2 
    64       REAL(wp) ::   zdispot, zfact, zcalcon, zalka, zaldi 
     63      INTEGER  ::   ji, jj, jk 
     64      REAL(wp) ::   zdispot, zfact, zcalcon 
    6565      REAL(wp) ::   zomegaca, zexcess, zexcess0 
    6666      CHARACTER (len=25) :: charout 
    67       REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss    
     67      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zco3sat, zcaldiss, zhinit, zhi    
    6868      !!--------------------------------------------------------------------- 
    6969      ! 
    7070      IF( nn_timing == 1 )  CALL timing_start('p4z_lys') 
    7171      ! 
    72       CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss ) 
     72      CALL wrk_alloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss, zhinit, zhi ) 
    7373      ! 
    7474      zco3    (:,:,:) = 0. 
    7575      zcaldiss(:,:,:) = 0. 
     76      zhinit(:,:,:)   = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn ) 
    7677      !     ------------------------------------------- 
    7778      !     COMPUTE [CO3--] and [H+] CONCENTRATIONS 
    7879      !     ------------------------------------------- 
    79        
    80       DO jn = 1, 5                               !  BEGIN OF ITERATION 
    81          ! 
    82 !CDIR NOVERRCHK 
    83          DO jk = 1, jpkm1 
    84 !CDIR NOVERRCHK 
    85             DO jj = 1, jpj 
    86 !CDIR NOVERRCHK 
    87                DO ji = 1, jpi 
    88                   zfact = rhop(ji,jj,jk) / 1000. + rtrn 
    89                   zph  = hi(ji,jj,jk) * tmask(ji,jj,jk) / zfact + ( 1.-tmask(ji,jj,jk) ) * 1.e-9 ! [H+] 
    90                   zdic  = trb(ji,jj,jk,jpdic) / zfact 
    91                   zalka = trb(ji,jj,jk,jptal) / zfact 
    92                   ! CALCULATE [ALK]([CO3--], [HCO3-]) 
    93                   zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph / ( aphscale(ji,jj,jk) + rtrn )  & 
    94                   &       + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 
    95                   ! CALCULATE [H+] and [CO3--] 
    96                   zaldi = zdic - zalk 
    97                   zah2  = SQRT( zaldi * zaldi + 4.* ( zalk * ak23(ji,jj,jk) / ak13(ji,jj,jk) ) * ( zdic + zaldi ) ) 
    98                   zah2  = 0.5 * ak13(ji,jj,jk) / zalk * ( zaldi + zah2 ) 
    99                   ! 
    100                   zco3(ji,jj,jk) = zalk / ( 2. + zah2 / ak23(ji,jj,jk) ) * zfact 
    101                   hi(ji,jj,jk)   = zah2 * zfact 
    102                END DO 
     80 
     81      CALL p4z_che_solve_hi( zhinit, zhi ) 
     82 
     83      DO jk = 1, jpkm1 
     84         DO jj = 1, jpj 
     85            DO ji = 1, jpi 
     86               zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   & 
     87               &                + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn ) 
     88               hi(ji,jj,jk)   = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000. 
    10389            END DO 
    10490         END DO 
    105          ! 
    106       END DO  
     91      END DO 
    10792 
    10893      !     --------------------------------------------------------- 
     
    138123              !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION 
    139124              zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution 
    140               zco3(ji,jj,jk)      = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk) 
    141125              ! 
    142126              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk) 
     
    167151      ENDIF 
    168152      ! 
    169       CALL wrk_dealloc( jpi, jpj, jpk, zco3, zco3sat, zcaldiss ) 
     153      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi, zco3sat ) 
    170154      ! 
    171155      IF( nn_timing == 1 )  CALL timing_stop('p4z_lys') 
     
    186170      !! 
    187171      !!---------------------------------------------------------------------- 
    188       INTEGER  ::  ji, jj, jk 
    189172      INTEGER  ::  ios                 ! Local integer output status for namelist read 
    190       REAL(wp) ::  zcaralk, zbicarb, zco3 
    191       REAL(wp) ::  ztmas, ztmas1 
    192  
    193173      NAMELIST/nampiscal/ kdca, nca 
    194174      !!---------------------------------------------------------------------- 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90

    r6420 r7607  
    8585        CALL p4z_che                              ! initialize the chemical constants 
    8686        ! 
    87         IF( .NOT. ln_rsttr ) THEN  ;   CALL p4z_ph_ini   !  set PH at kt=nit000  
     87        IF( .NOT. ln_rsttr ) THEN  ;   CALL p4z_che_ahini( hi )   !  set PH at kt=nit000 
    8888        ELSE                       ;   CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields  
    8989        ENDIF 
     
    310310   END SUBROUTINE p4z_sms_init 
    311311 
    312    SUBROUTINE p4z_ph_ini 
    313       !!--------------------------------------------------------------------- 
    314       !!                   ***  ROUTINE p4z_ini_ph  *** 
    315       !! 
    316       !!  ** Purpose : Initialization of chemical variables of the carbon cycle 
    317       !!--------------------------------------------------------------------- 
    318       INTEGER  ::  ji, jj, jk 
    319       REAL(wp) ::  zcaralk, zbicarb, zco3 
    320       REAL(wp) ::  ztmas, ztmas1 
    321       !!--------------------------------------------------------------------- 
    322  
    323       ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???) 
    324       ! -------------------------------------------------------- 
    325       DO jk = 1, jpk 
    326          DO jj = 1, jpj 
    327             DO ji = 1, jpi 
    328                ztmas   = tmask(ji,jj,jk) 
    329                ztmas1  = 1. - tmask(ji,jj,jk) 
    330                zcaralk = trb(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  ) 
    331                zco3    = ( zcaralk - trb(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1 
    332                zbicarb = ( 2. * trb(ji,jj,jk,jpdic) - zcaralk ) 
    333                hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1 
    334             END DO 
    335          END DO 
    336      END DO 
    337      ! 
    338    END SUBROUTINE p4z_ph_ini 
    339  
    340312   SUBROUTINE p4z_rst( kt, cdrw ) 
    341313      !!--------------------------------------------------------------------- 
     
    350322      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
    351323      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    352       ! 
    353       INTEGER  ::  ji, jj, jk 
    354       REAL(wp) ::  zcaralk, zbicarb, zco3 
    355       REAL(wp) ::  ztmas, ztmas1 
    356324      !!--------------------------------------------------------------------- 
    357325 
     
    365333            CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  ) 
    366334         ELSE 
    367 !            hi(:,:,:) = 1.e-9  
    368             CALL p4z_ph_ini 
     335            CALL p4z_che_ahini( hi ) 
    369336         ENDIF 
    370337         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) ) 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90

    r6287 r7607  
    9393 
    9494   !!* Variable for chemistry of the CO2 cycle 
    95    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ??? 
    9695   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak13       !: ??? 
    9796   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak23       !: ??? 
    9897   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksp       !: ??? 
    99    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akw3       !: ??? 
    100    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   borat      !: ??? 
    10198   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi         !: ??? 
    10299   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   excess     !: ??? 
     
    153150 
    154151      !* Variable for chemistry of the CO2 cycle 
    155       ALLOCATE( akb3(jpi,jpj,jpk)    , ak13  (jpi,jpj,jpk) ,       & 
     152      ALLOCATE( ak13  (jpi,jpj,jpk) ,                              & 
    156153         &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,       & 
    157          &      akw3(jpi,jpj,jpk)    , borat (jpi,jpj,jpk) ,       & 
    158154         &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,       & 
    159155         &      aphscale(jpi,jpj,jpk),                           STAT=ierr(4) ) 
  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90

    r7522 r7607  
    131131      INTEGER, INTENT(in) ::   kt 
    132132      INTEGER  :: jn 
    133       REAL(wp) :: zkt 
     133      REAL(wp) :: zkt, zrec 
    134134      CHARACTER(len=1)               ::   cl1                      ! 1 character 
    135135      CHARACTER(len=2)               ::   cl2                      ! 2 characters 
     
    153153         ! 
    154154         !                                            !* Restart: read in restart file 
    155          IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 .AND. & 
    156                             iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 .AND. & 
    157                             iom_varid( numrtr, 'ktdcy'    , ldstop = .FALSE. ) > 0 ) THEN  
    158             CALL iom_get( numrtr, 'ktdcy', zkt )   !  A mean of qsr 
     155         IF(  ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0  & 
     156           &           .AND. iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0  & 
     157           &           .AND. iom_varid( numrtr, 'ktdcy'    , ldstop = .FALSE. ) > 0  & 
     158           &           .AND. iom_varid( numrtr, 'nrdcy'    , ldstop = .FALSE. ) > 0  ) THEN 
     159 
     160            CALL iom_get( numrtr, 'ktdcy', zkt )   
    159161            rsecfst = INT( zkt ) * rdttrc(1) 
    160162            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s ' 
    161163            CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean )   !  A mean of qsr 
    162             DO jn = 1, nb_rec_per_day  
    163              IF( jn <= 9 )  THEN 
    164                WRITE(cl1,'(i1)') jn 
    165                CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr 
    166              ELSE 
    167                WRITE(cl2,'(i2.2)') jn 
    168                CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr 
    169              ENDIF 
    170            ENDDO 
     164            CALL iom_get( numrtr, 'nrdcy', zrec )   !  Number of record per days 
     165            IF( INT( zrec ) == nb_rec_per_day ) THEN 
     166               DO jn = 1, nb_rec_per_day  
     167                  IF( jn <= 9 )  THEN 
     168                    WRITE(cl1,'(i1)') jn 
     169                    CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr 
     170                  ELSE 
     171                    WRITE(cl2,'(i2.2)') jn 
     172                    CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr 
     173                  ENDIF 
     174              ENDDO 
     175            ELSE 
     176               DO jn = 1, nb_rec_per_day 
     177                  qsr_arr(:,:,jn) = qsr_mean(:,:) 
     178               ENDDO 
     179            ENDIF 
    171180         ELSE                                         !* no restart: set from nit000 values 
    172181            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean set to nit000 values' 
     
    185194      llnew   = ( rseclast - rsecfst ) .ge.  rdt_sampl    !   new shortwave to store 
    186195      IF( llnew ) THEN 
    187           IF( lwp ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, & 
     196          IF( lwp .AND. kt < nittrc000 + 100 ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, & 
    188197             &                      ' time = ', rseclast/3600.,'hours ' 
    189198          rsecfst = rseclast 
     
    199208         IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file  kt =', kt 
    200209         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    201          zkt = REAL( kt, wp ) 
    202          CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt ) 
     210         zkt  = REAL( kt, wp ) 
     211         zrec = REAL( nb_rec_per_day, wp ) 
     212         CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt  ) 
     213         CALL iom_rstput( kt, nitrst, numrtw, 'nrdcy', zrec ) 
    203214          DO jn = 1, nb_rec_per_day  
    204215             IF( jn <= 9 )  THEN 
Note: See TracChangeset for help on using the changeset viewer.