Changeset 27 for branches/2016


Ignore:
Timestamp:
10/06/16 11:31:57 (8 years ago)
Author:
vancop
Message:

First physical fixes, new ice liquid fraction and mushy-layer thermal conductivity

Location:
branches/2016/dev_v3.20_2016_platelet
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • branches/2016/dev_v3.20_2016_platelet/INPUT/ISPOL/source_3.20/icephys.param

    r6 r27  
    5151# delta_cw      ! CW88 grav. drainage eff.(1.68e-7 for CW88 or 0.588e-7 for V08 value) 
    52521.68e-7 
     53# rn_e_newice   ! Liquid fraction in new ice (max fraction due to flottability of ice) 
     540.75 
     55# nn_thcond     ! Type of thermal conductivity (1- Pringle et al 2007; 2- Mushy layer) 
     562 
    5357# 
    5458#-------------------------------------------------------------------------------------------------- 

  • branches/2016/dev_v3.20_2016_platelet/SCRIPTS/lim1d_ispol.nqs

    r6 r27  
    99USERNICK="Martouf"            # user nickname 
    1010 
    11 EXP_ID="test_3.20_ISPOL_a"      # name of the experiment ### CHANGE 
     11EXP_ID="test_3.20_ISPOL_01"   # name of the experiment ### CHANGE 
    1212SOURCE="source_3.20"          # version of the sources 
    1313TS="1d"                       # time step (1h or 1d) 
    … …  
    1717SUBDIR=`pwd` 
    1818 
    19 MAINDIR="$HOME/Boulot/CODES/LIM1D/2015_v3.20" ### CHANGE 
     19MAINDIR="$HOME/Boulot/CODES/LIM1D_SVN/dev_v3.20_2016_platelet" ### CHANGE 
    2020 
    2121SCRATCHDIR="$MAINDIR/SCRATCH/$EXP_ID"        # SCRATCH - temporary directory on which code is run 
    … …  
    2424INFILEDIR="$MAINDIR/INPUT/$CONF/$SOURCE"     # INPUT - initialization, forcing & param files 
    2525 
    26 GRAPHDIR="$HOME/Boulot/SCIENCE/PLOT_SCRIPTS/LIM1D_BIO/IDL/$CONF" # IDL plots 
     26GRAPHDIR="$HOME/Boulot/SCIENCE/PLOT_SCRIPTS/LIM1D/IDL/$CONF" # IDL plots 
    2727 
    2828####################################################################################### 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice.com

    r6 r27  
    181181 
    182182      COMMON/barrowconf/sal_read(11), hi_read(11), hgins, hnins, 
    183      &              tsuins, oce_sal, oce_flx, num_sal, nday1, 
     183     &              tsuins, s_w, oce_flx, num_sal, nday1, 
    184184     &              i_sal 
    185185      
    … …  
    191191      COMMON/fluidtpt/flu_beta, rad_io, 
    192192     &                frtr_si_phy, qsummer, d_br_mol, d_br_tur, 
    193      &                ra_c, ra_smooth, e_thr_flu, delta_cw, ini_swi,  
     193     &                ra_c, ra_smooth, e_thr_flu, delta_cw,  
     194     &                rn_e_newice, ini_swi,  
     195     &                nn_thcond, 
    194196     &                s_ini, ln_flo, ln_flu, ln_grd 
    195197 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_bio_remap.f

    r6 r27  
    169169 
    170170         IF ( ln_trbo )  
    171      &      ch_bo_bio(jn) = e_skel * c_w_bio(jn) * MAX( dh_i_bott(ji), 
    172      &                      0.0 ) 
     171     &      ch_bo_bio(jn) = rn_e_newice * c_w_bio(jn) * 
     172     &                      MAX( dh_i_bott(ji), 0.0 ) 
    173173 
    174174         IF ( ln_trsi ) 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_brine.f

    r6 r27  
    4848         sbr_bio(layer)   = -21.4*tc_bio(layer) - 0.886*zt2 -0.0107*zt3 ! Brine 3rd order salinity 
    4949          
    50          rhobr_bio(layer) = rho0 + beta_ocs*( sbr_bio(layer)-oce_sal )  ! Brine density 
     50         rhobr_bio(layer) = rho0 + beta_ocs*( sbr_bio(layer)-s_w )  ! Brine density 
    5151             
    5252      END DO 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_output.f

    r6 r27  
    208208         WRITE(numout,*) ' fsbp dimension created ' 
    209209 
    210          CALL CF_CREATE_VAR("oce_sal","Ocean Salinity", 
     210         CALL CF_CREATE_VAR("s_w","Ocean Salinity", 
    211211     .                   "kgNaCl/dm3","time","-","-","-" ) 
    212          WRITE(numout,*) ' oce_sal dimension created ' 
     212         WRITE(numout,*) ' s_w dimension created ' 
    213213 
    214214         ! Forcing 
    … …  
    590590     &               1, 1, 1, REAL(fsbp) ) 
    591591 
    592       CALL CF_WRITE (filenc, 'oce_sal', numit-nstart+1,  
    593      &               1, 1, 1, REAL(oce_sal) ) 
     592      CALL CF_WRITE (filenc, 's_w', numit-nstart+1,  
     593     &               1, 1, 1, REAL(s_w) ) 
    594594 
    595595      ! Forcing 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_phy_ini.f

    r6 r27  
    8686      CALL CF_CLOSE (filenc) ! close nc file 
    8787 
    88       t_bo_b(ji)   =  273.15-tmut*oce_sal  
     88      t_bo_b(ji)   =  273.15-tmut*s_w 
    8989 
    9090      !---------------------- 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_phy_param.f

    r6 r27  
    108108      rho0    = 1025.0            ! ocean mean density 
    109109      cpw     = 3.99d+03          ! seawater specific heat 
    110       oce_sal = 34.0              ! 34. is the control 
     110      s_w    = 34.0              ! 34. is the control 
    111111 
    112112      visc_br = 1.79e-3           ! dynamic viscosity of water at 0C  
    … …  
    191191      READ(25,*) delta_cw         ! Cox and weeks gravity drainage parameter 
    192192      READ(25,*) 
     193      READ(25,*) rn_e_newice      ! Liquid fraction in new ice 
     194      READ(25,*) 
     195      READ(25,*) nn_thcond        ! Type of thermal conductivity computation 
     196      READ(25,*) 
    193197      READ(25,*) 
    194198      READ(25,*) 
    … …  
    251255      WRITE(numout,*) ' gravdr     :', gravdr 
    252256      WRITE(numout,*) ' delta_cw   :', delta_cw 
     257      WRITE(numout,*) ' rn_e_newice:', rn_e_newice 
     258      WRITE(numout,*) ' nn_thcond  :', nn_thcond 
    253259      WRITE(numout,*) ' emig        :', emig 
    254260      WRITE(numout,*) ' fpar_fsw    :', fpar_fsw 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_phy_remap.f

    r6 r27  
    436436! 
    437437      ! basal new ice salt content 
    438       zsh_i_new = e_skel * oce_sal * MAX ( dh_i_bott(ji), 0.0 ) 
     438      zsh_i_new = rn_e_newice * s_w * MAX ( dh_i_bott(ji), 0.0 ) 
    439439      ! snow ice salt content 
    440       s_i_snic  =  ( rhog - rhon ) / rhog * oce_sal *  
     440      s_i_snic  =  ( rhog - rhon ) / rhog * s_w * 
    441441     &              frtr_si_phy ! snow ice salinity 
    442442      zsh_i_sni = s_i_snic * dh_snowice(ji) 
    … …  
    470470      IF ( ln_write ) THEN 
    471471         WRITE(numout,*) ' frtr_si_phy:', frtr_si_phy 
    472          WRITE(numout,*) ' oce_sal  : ', oce_sal 
     472         WRITE(numout,*) ' s_w      : ', s_w 
    473473         WRITE(numout,*) ' s_i_snic : ', s_i_snic 
    474474         WRITE(numout,*) ' zsh_i0   : ', zsh_i0(ntop0:nbot0) 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_sal_diff.f

    r6 r27  
    194194      !-------------------- 
    195195      DO layer = 1, nlay_i   
    196          zrhodiff(layer) = - beta_ocs * ( oce_sal - zsigma(layer) ) 
     196         zrhodiff(layer) = - beta_ocs * ( s_w - zsigma(layer) ) 
    197197      END DO 
    198198 
    … …  
    444444      ztrid(nlay_i,3) = 0. 
    445445      zind(nlay_i)    = z_sbr_i(nlay_i) + za(nlay_i) * ( zb(nlay_i)  
    446      &                - zc(nlay_i) ) * oce_sal 
     446     &                - zc(nlay_i) ) * s_w 
    447447 
    448448      IF ( ln_write ) THEN 
    … …  
    523523      zfb      =  - e_i_b(nlay_i) *  
    524524     &             ( diff_br(nlay_i) * 2.0 / deltaz_i_phy(nlay_i) *  
    525      &             ( z_sbr_i(nlay_i) - oce_sal ) + w_flood * ( z_flood * 
    526      &             oce_sal + ( 1. - z_flood ) * z_sbr_i(nlay_i) ) ) 
     525     &             ( z_sbr_i(nlay_i) - s_w ) + w_flood * ( z_flood * 
     526     &             s_w + ( 1. - z_flood ) * z_sbr_i(nlay_i) ) ) 
    527527     &            - e_i_b(nlay_i) * w_flush * z_sbr_i(nlay_i) 
    528528!    &            - qsummer * z_sbr_i(nlay_i) / ddtb 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_sal_diff_CW.f

    r6 r27  
    429429     &                * diff_br(nlay_i) * 2.0 
    430430     &                / deltaz_i_phy(nlay_i) * ( z_sbr_i(nlay_i) 
    431      &                  - oce_sal ) ) * zswitch_open 
     431     &                  - s_w ) ) * zswitch_open 
    432432     &                + zswitchs * ( - qsummer * z_sbr_i(nlay_i) ) 
    433433     &                / ddtb 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_th_dh.f

    r6 r27  
    315315      ! Basal growth 
    316316      !-------------- 
    317       ! Initial value (tested 1D, can be anything between 1 and 20) 
    318       s_i_new      = 10.0 
    319       num_iter_max = 10 
    320  
    321       ! the growth rate (dh_i_bott) is function of the new ice 
    322       ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice 
    323       ! salinity (snewice). snewice depends on dh_i_bott 
    324       ! it converges quickly, so, no problem 
    325  
    326       ! Iterative procedure 
    327317      IF ( ( fc_bo_i(ji) + fbbqb(ji) ) .LT. 0.0 ) THEN 
    328  
    329          IF ( l_write ) WRITE(numout,*) ' Energy available for  
    330      &                                    basal growth : ', 
    331      &                  fc_bo_i(ji) + fbbqb(ji) 
    332  
    333          DO iter = 1, num_iter_max 
    334                tmelts             =   - tmut * s_i_new + tpw ! Melting point in K 
    335                q_i_b(ji,nlay_i+1) = rhog*( cpg*(tmelts-t_bo_b(ji)) 
    336      &                 + lfus*( 1.0-(tmelts-tpw) 
    337      &                 / (t_bo_b(ji) - tpw) ) 
    338      &                 - cpw*(tmelts-tpw) )  
    339  
    340                ! Bottom growth rate = - F*dt / q 
    341                dh_i_bott(ji)    =  - ddtb*(fc_bo_i(ji) + fbbqb(ji) ) 
    342      &                             / q_i_b(ji,nlay_i+1) 
    343                ! 
    344                ! New ice salinity ( Cox and Weeks, JGR, 1988 ) 
    345                ! 
    346                ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7 
    347                ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 
    348                ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8 
    349                ! 
    350                zgrr   = MIN( 1.0e-3,  
    351      &                  MAX ( dh_i_bott(ji) / ddtb , zeps ) ) 
    352                zswi2  = MAX( 0.0 , SIGN( 1.0d0 , zgrr - 3.6e-7 ) ) 
    353                zswi12 = MAX( 0.0 , SIGN( 1.0d0 , zgrr - 2.0e-8 ) ) *  
    354      &                  ( 1.0 - zswi2 ) 
    355                zswi1  = 1. - zswi2 * zswi12 
    356                zfracs = zswi1  * 0.12 +    
    357      &                  zswi12 * ( 0.8925 + 0.0568 *  
    358      &                  LOG( 100.0 * zgrr ) ) +   
    359      &                  zswi2  * 0.26 /    
    360      &                  ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
    361                zfracs = 1. 
    362                zds     = zfracs * oce_sal - s_i_new 
    363                s_i_new = zfracs * oce_sal 
    364  
    365                ! the brine volume in the skeletal layer is equal to f 
    366                e_skel  = zfracs 
    367  
    368                ! salt flux due to initial salt entrapment 
    369                fsbp = oce_sal * ( 1. - zfracs ) * dh_i_bott(ji) / ddtb * 
    370      &                rhog / 1000. 
    371  
    372 !              WRITE(numout,*) 
    373 !              WRITE(numout,*) ' ddtb      : ', ddtb 
    374 !              WRITE(numout,*) ' zgrr      : ', zgrr 
    375 !              WRITE(numout,*) ' zswi12    : ', zswi12 
    376 !              WRITE(numout,*) ' zswi1     : ', zswi1 
    377 !              WRITE(numout,*) ' zswi2     : ', zswi2 
    378 !              WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 
    379 !              WRITE(numout,*) ' oce_sal   : ', oce_sal 
    380 !              WRITE(numout,*) ' zfracs    : ', zfracs  
    381 !              WRITE(numout,*) ' s_i_new   : ', s_i_new 
    382 !              WRITE(numout,*) 
    383  
    384          END DO ! iter 
    385       ENDIF ! fc_bo_i 
    386  
    387       IF ( l_write ) THEN 
    388          WRITE(numout,*) ' e_skel : ', e_skel 
    389          WRITE(numout,*) 
    390       ENDIF 
    391  
    392       ! Final values 
    393       IF ( (fc_bo_i(ji)+fbbqb(ji)) .LT. 0.0 ) THEN 
    394       ! New ice salinity must not exceed 15 psu 
    395          zoldsinew   = s_i_new 
    396 ! TEST MARTIN NOV 2010 this line should be removed in future versions 
    397          s_i_new     = MIN( s_i_new , s_i_max ) 
    398          ! Metling point in K 
    399          tmelts             =   - tmut * s_i_new + tpw 
     318  
     319         s_i_new = rn_e_newice * s_w ! New ice salinity (g/kg) 
     320  
     321         ztmelts  = - tmut * s_i_new + tpw ! Melting point in K 
     322   
    400323         ! New ice heat content (Bitz and Lipscomb, 1999) 
     324         ! should be zdE instead 
    401325         q_i_b(ji,nlay_i+1) = rhog*( cpg*(tmelts-t_bo_b(ji)) 
    402326     &                 + lfus*( 1.0-(tmelts-tpw) 
    403327     &                 / (t_bo_b(ji) - tpw) ) 
    404328     &                 - cpw*(tmelts-tpw) )  
     329 
     330         zE1 = - cpw * t_bo_b(ji) ! specific enthalpy of sea water 
     331         zE2 = q_i_b(ji,nlay_i+1) ! specific enthalpy of new sea ice 
     332 
     333         zdE = zE2 - zE1 
     334 
    405335         dh_i_bott(ji)    =  - ddtb*(fc_bo_i(ji) + fbbqb(ji) ) 
    406      &                    / q_i_b(ji,nlay_i+1) 
     336     &                    / zdE 
     337   
     338         ! salt flux due to initial salt entrapment (keep ?) 
     339         fsbp = s_w * ( 1. - rn_e_newice ) * dh_i_bott(ji) / ddtb * 
     340     &          rhog / 1000. 
     341   
    407342         IF ( l_write ) WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 
    408  
    409       ENDIF  
     343   
     344      ENDIF 
    410345 
    411346      !----------------- 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/ice_th_diff.f

    r6 r27  
    5959      INTEGER    numeqmin, numeqmax, numeq 
    6060      DIMENSION  ztcond_i(0:nlay_i), 
     61     &           ze(0:nlay_i), 
     62     &           ztc(0:nlay_i), 
    6163     &           zkappa_s(0:nlay_s), 
    6264     &           zkappa_i(0:nlay_i),ztstemp(0:nlay_s),ztitemp(0:nlay_i), 
    … …  
    122124!------------------------------------------------------------------------------  
    123125! 
    124       ! Pringle et al., JGR 2007 formula 
    125       ! 2.11 + 0.09 S/T - 0.011.T 
     126      nn_thcond = 2 
    126127 
    127128      DO 10 ji = kideb, kiut 
    128129 
    129       ! thermal conductivity in the snow 
     130      ! thermal conductivity of the snow 
    130131      ykn(ji)  =  xkn 
    131132      zkimin   =  0.1 
     133 
     134      ! thermal conductivity of the sea ice 
     135 
     136      SELECT CASE ( nn_thcond ) 
     137 
     138      !--------------------------------------------------- 
     139      CASE(1) ! Pringle et al., JGR 2007 
     140              ! 2.11 + 0.09 S/T - 0.011.T 
     141              ! Empirical, direct measurements in sea ice 
     142      !--------------------------------------------------- 
    132143 
    133144      ztcond_i(0)     = xkg + betak1*s_i_b(ji,1)  
    134145     &                   / MIN( -zeps , t_i_b(ji,1) - tpw )  
    135146     &                   - betak2* ( t_i_b(ji,1) - tpw ) 
    136       ztcond_i(0)     = MAX( ztcond_i(0) , zkimin ) 
    137147 
    138148      DO layer = 1, nlay_i-1 
    … …  
    142152     &                      - betak2 * 0.5 * ( t_i_b(ji,layer) +  ! bugfix fred dupont add 0.5 
    143153     &                        t_i_b(ji,layer+1) - 2.0*tpw ) 
    144          ztcond_i(layer) = MAX( ztcond_i(layer) , zkimin ) 
    145154      END DO 
    146155 
    … …  
    148157     &                        MIN( -zeps , t_bo_b(ji) - tpw )  
    149158     &                      - betak2 * ( t_bo_b(ji) - tpw ) 
     159 
     160      !--------------------------------------------------- 
     161      CASE(2) ! Mushy layer: k = ki . ( 1 - phi) + kbr . phi 
     162              ! Thermal conductivity of pure ice at interfaces from Yen (1991) 
     163              ! zki = 2.21 - 1.0e-2*ztc + 3.44e-5*ztc.^2 (Yen 1991)  
     164              ! Thermal conductivity of brine at interfaces (Castelli et al DSR 1974) 
     165              ! 0.55286 + 1.8364e-3 * TT - 3.3058e-7 * TT.^3 
     166      !--------------------------------------------------- 
     167 
     168      ! Brine fraction at layer mid-points 
     169      e_i_b(:) = -tmut * s_i_b(ji,:) / ( t_i_b(ji,:) - tpw ) 
     170 
     171      ! Brine fraction at interfaces 
     172      zde_dz       = ( e_i_b(2)   - e_i_b(1) ) /  
     173     &               ( z_i_phy(2) - z_i_phy(1)  ) 
     174      ze(0)        = e_i_b(1)     - zde_dz * deltaz_i_phy(1) / 2. 
     175 
     176      DO layer = 1, nlay_i - 1 
     177         zde_dz    = ( e_i_b(layer+1) - e_i_b(layer) ) /  
     178     &               ( z_i_phy(layer+1) - z_i_phy(layer)  ) 
     179         ze(layer) = e_i_b(layer)    + zde_dz * deltaz_i_phy(layer) / 2. 
     180      END DO 
     181 
     182      ze(nlay_i) = 1. 
     183 
     184      ! Temperature at interfaces 
     185      !!! ztc(0) could be better interpolated ... use effective conductivity ? 
     186      zdt_dz       = ( t_i_b(ji,2)     - t_i_b(ji,1) ) /  
     187     &               ( z_i_phy(2)      - z_i_phy(1)  ) 
     188      ztc(0)       = t_i_b(ji,1)       - zdt_dz * deltaz_i_phy(1) / 2. 
     189      ztc(0)       = ztc(0) - tpw 
     190 
     191      DO layer = 1, nlay_i - 1 
     192         zdt_dz    = ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) /  
     193     &               ( z_i_phy(layer+1)  - z_i_phy(layer)  ) 
     194         ztc(layer) = t_i_b(ji,layer) +  
     195     &                zdt_dz * deltaz_i_phy(layer) / 2. 
     196         ztc(layer) = ztc(layer) - tpw 
     197      END DO 
     198      ztc(nlay_i)    = t_bo_b(ji) - tpw 
     199 
     200      DO layer = 0, nlay_i 
     201         zt1 = ztc(layer) 
     202         zt2 = zt1*zt1 
     203         zt3 = zt1*zt2 
     204         zki = 2.2156 - 1.0045e-2*zt1 + 3.44520e-5*ztc2 
     205         zkbr = 0.55286 + 1.8364e-3*zt1 - 3.3058e-7*zt3 
     206         ztcond_i(layer) = zki * ( 1. - ze(layer) ) + zkbr * ze(layer) 
     207      END DO 
     208 
     209      END SELECT 
     210 
     211      ! Cap thermal conductivity to prevent small values 
     212      ztcond_i(0)     = MAX( ztcond_i(0) , zkimin ) 
     213      DO layer = 0, nlay_i 
     214         ztcond_i(layer) = MAX( ztcond_i(layer) , zkimin ) 
     215      END DO 
    150216      ztcond_i(nlay_i)   = MAX( ztcond_i(nlay_i) , zkimin ) 
    151217 
    … …  
    488554         zindterm(numeqmin) =  ztiold(1) + zeta_i(1)* 
    489555!    &                         (radab_phy_i(1) + radab_alg_i(1) 
    490      &                         ( radab_i(1) +  
     556     &                         ( radab_i(1)    
    491557     &                         + zkappa_i(1)*t_bo_b(ji) )  
    492558     &                         + t_su_b(ji)*zeta_i(1)*zkappa_i(0)*2.0 

  • branches/2016/dev_v3.20_2016_platelet/SOURCES/source_3.20/thermo.com

    r6 r27  
    9999     &                 q_i_layer_fin(nbpt,maxnlay), fprec, fsnic 
    100100 
    101       common/ salt / beta_sal, s_i_new, s_i_snic, e_skel, q_summer, 
     101      common/ salt / beta_sal, s_i_new, s_i_snic, q_summer, 
    102102     &               diff_br(maxnlay), rayleigh(maxnlay), fsb, fsbp, 
    103103     &               w_flood, w_flush 
Note: See TracChangeset for help on using the changeset viewer.