- Timestamp:
- 2019-08-02T16:19:00+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/ENHANCE-02_ISF_nemo
- Files:
-
- 1 added
- 2 copied
-
. (copied) (copied from NEMO/trunk)
-
src/OCE/ISF (added)
-
src/OCE/ISF/isfmlt.F90 (copied) (copied from NEMO/trunk/src/OCE/SBC/sbcisf.F90) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfmlt.F90
r11310 r11395 1 MODULE sbcisf1 MODULE isfmlt 2 2 !!====================================================================== 3 3 !! *** MODULE sbcisf *** 4 !! Surface module : update surface ocean boundary condition under ice 5 !! shelf 4 !! Surface module : compute iceshelf melt and heat flux 6 5 !!====================================================================== 7 6 !! History : 3.2 ! 2011-02 (C.Harris ) Original code isf cav 8 7 !! X.X ! 2006-02 (C. Wang ) Original code bg03 9 8 !! 3.4 ! 2013-03 (P. Mathiot) Merging + parametrization 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! sbc_isf : update sbc under ice shelf 9 !! 4.1 ! 2019-09 (P. Mathiot) Split param/explicit ice shelf and re-organisation 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! isfmlt : compute iceshelf melt and heat flux 14 14 !!---------------------------------------------------------------------- 15 15 USE oce ! ocean dynamics and tracers … … 25 25 USE lbclnk ! 26 26 USE lib_fortran ! glob_sum 27 ! 28 USE isfrst ! iceshelf restart 29 USE isftbl ! ice shelf boundary layer 30 USE isfpar ! ice shelf parametrisation 31 USE isfcav ! ice shelf cavity 32 USE isf ! isf variables 27 33 28 34 IMPLICIT NONE 35 29 36 PRIVATE 30 37 31 PUBLIC sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc ! routine called in sbcmod and divhor 32 33 ! public in order to be able to output then 34 35 REAL(wp), PUBLIC :: rn_hisf_tbl !: thickness of top boundary layer [m] 36 INTEGER , PUBLIC :: nn_isf !: flag to choose between explicit/param/specified 37 INTEGER , PUBLIC :: nn_isfblk !: flag to choose the bulk formulation to compute the ice shelf melting 38 INTEGER , PUBLIC :: nn_gammablk !: flag to choose how the exchange coefficient is computed 39 REAL(wp), PUBLIC :: rn_gammat0 !: temperature exchange coeficient [] 40 REAL(wp), PUBLIC :: rn_gammas0 !: salinity exchange coeficient [] 41 42 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt , misfkb !: Level of ice shelf base 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rzisf_tbl !: depth of calving front (shallowest point) nn_isf ==2/3 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl, rhisf_tbl_0 !: thickness of tbl [m] 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hisf_tbl !: 1/thickness of tbl 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ralpha !: proportion of bottom cell influenced by tbl 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfLeff !: effective length (Leff) BG03 nn_isf==2 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ttbl, stbl, utbl, vtbl !: top boundary layer variable at T point 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf [W/m2] 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc !: before and now T & S isf contents [K.m/s & PSU.m/s] 51 52 LOGICAL, PUBLIC :: l_isfcpl = .false. !: isf recieved from oasis 53 54 REAL(wp), PUBLIC, SAVE :: rcpisf = 2000.0_wp !: specific heat of ice shelf [J/kg/K] 55 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp !: heat diffusivity through the ice-shelf [m2/s] 56 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp !: volumic mass of ice shelf [kg/m3] 57 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp !: air temperature on top of ice shelf [C] 58 REAL(wp), PUBLIC, SAVE :: rLfusisf = 0.334e6_wp !: latent heat of fusion of ice shelf [J/kg] 59 60 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 61 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 62 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf melting file to be read 63 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 64 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the isf melting param. file to be read 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 66 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf !: information about the grounding line depth file to be read 67 TYPE(FLD_N) , PUBLIC :: sn_depmin_isf !: information about the calving line depth file to be read 68 TYPE(FLD_N) , PUBLIC :: sn_Leff_isf !: information about the effective length file to be read 69 38 PUBLIC isf_mlt, isf_mlt_init ! routine called in sbcmod and divhor 39 70 40 !!---------------------------------------------------------------------- 71 41 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 75 45 CONTAINS 76 46 77 SUBROUTINE sbc_isf( kt )47 SUBROUTINE isf_mlt( kt ) 78 48 !!--------------------------------------------------------------------- 79 !! *** ROUTINE sbc_isf *** 80 !! 81 !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 82 !! melting and freezing 83 !! 84 !! ** Method : 4 parameterizations are available according to nn_isf 85 !! nn_isf = 1 : Realistic ice_shelf formulation 86 !! 2 : Beckmann & Goose parameterization 87 !! 3 : Specified runoff in deptht (Mathiot & al. ) 88 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 49 !! *** ROUTINE isf_mlt *** 50 !! 51 !! ** Purpose : 52 !! 53 !! ** Method : 54 !! 89 55 !!---------------------------------------------------------------------- 90 56 INTEGER, INTENT(in) :: kt ! ocean time step 91 !92 INTEGER :: ji, jj, jk ! loop index93 INTEGER :: ikt, ikb ! local integers94 REAL(wp), DIMENSION(jpi,jpj) :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)95 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zqhcisf2d96 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfwfisf3d, zqhcisf3d, zqlatisf3d97 57 !!--------------------------------------------------------------------- 98 58 ! 99 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN ! compute salt and heat flux 100 ! 101 SELECT CASE ( nn_isf ) 102 CASE ( 1 ) ! realistic ice shelf formulation 103 ! compute T/S/U/V for the top boundary layer 104 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 105 CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 106 CALL sbc_isf_tbl(un(:,:,:) ,utbl(:,:),'U') 107 CALL sbc_isf_tbl(vn(:,:,:) ,vtbl(:,:),'V') 108 ! iom print 109 CALL iom_put('ttbl',ttbl(:,:)) 110 CALL iom_put('stbl',stbl(:,:)) 111 CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 112 CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)) 113 ! compute fwf and heat flux 114 ! compute fwf and heat flux 115 IF( .NOT.l_isfcpl ) THEN ; CALL sbc_isf_cav (kt) 116 ELSE ; qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux 117 ENDIF 118 ! 119 CASE ( 2 ) ! Beckmann and Goosse parametrisation 120 stbl(:,:) = soce 121 CALL sbc_isf_bg03(kt) 122 ! 123 CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation) 124 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 125 IF( .NOT.l_isfcpl ) THEN 126 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 127 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fresh water flux from the isf (fwfisf <0 mean melting) 128 ENDIF 129 qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux 130 stbl(:,:) = soce 131 ! 132 CASE ( 4 ) ! specified fwf and heat flux forcing beneath the ice shelf 133 ! ! specified fwf and heat flux forcing beneath the ice shelf 134 IF( .NOT.l_isfcpl ) THEN 135 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 136 !CALL fld_read ( kt, nn_fsbc, sf_qisf ) 137 fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1) ! fwf 138 ENDIF 139 qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux 140 stbl(:,:) = soce 141 ! 142 END SELECT 143 144 ! compute tsc due to isf 145 ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 146 ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0). 147 ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 148 DO jj = 1,jpj 149 DO ji = 1,jpi 150 zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj)) 151 END DO 152 END DO 153 CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 154 155 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 156 risf_tsc(:,:,jp_sal) = 0.0_wp 157 158 ! lbclnk 159 CALL lbc_lnk_multi( 'sbcisf', risf_tsc(:,:,jp_tem), 'T', 1., risf_tsc(:,:,jp_sal), 'T', 1., fwfisf,'T', 1., qisf, 'T', 1.) 160 ! output 161 IF( iom_use('iceshelf_cea') ) CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) ) ! isf mass flux 162 IF( iom_use('hflx_isf_cea') ) CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rau0 * rcp ) ! isf sensible+latent heat (W/m2) 163 IF( iom_use('qlatisf' ) ) CALL iom_put( 'qlatisf' , qisf(:,:) ) ! isf latent heat 164 IF( iom_use('fwfisf' ) ) CALL iom_put( 'fwfisf' , fwfisf(:,:) ) ! isf mass flux (opposite sign) 165 166 ! Diagnostics 167 IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 168 ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) ) 169 ALLOCATE( zqhcisf2d(jpi,jpj) ) 170 ! 171 zfwfisf3d (:,:,:) = 0._wp ! 3d ice shelf melting (kg/m2/s) 172 zqhcisf3d (:,:,:) = 0._wp ! 3d heat content flux (W/m2) 173 zqlatisf3d(:,:,:) = 0._wp ! 3d ice shelf melting latent heat flux (W/m2) 174 zqhcisf2d (:,:) = fwfisf(:,:) * zt_frz * rcp ! 2d heat content flux (W/m2) 175 ! 176 DO jj = 1,jpj 177 DO ji = 1,jpi 178 ikt = misfkt(ji,jj) 179 ikb = misfkb(ji,jj) 180 DO jk = ikt, ikb - 1 181 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 182 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 183 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 184 END DO 185 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) & 186 & * ralpha(ji,jj) * e3t_n(ji,jj,jk) 187 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) & 188 & * ralpha(ji,jj) * e3t_n(ji,jj,jk) 189 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) & 190 & * ralpha(ji,jj) * e3t_n(ji,jj,jk) 191 END DO 192 END DO 193 ! 194 CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:)) 195 CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:)) 196 CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:)) 197 CALL iom_put('qhcisf' , zqhcisf2d (:,: )) 198 ! 199 DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d ) 200 DEALLOCATE( zqhcisf2d ) 201 ENDIF 202 ! 203 ENDIF 204 205 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 206 IF( ln_rstart .AND. & ! Restart: read in restart file 207 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 208 IF(lwp) WRITE(numout,*) ' nit000-1 isf tracer content forcing fields read in the restart file' 209 CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) , ldxios = lrxios ) ! before salt content isf_tsc trend 210 CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b' , risf_tsc_b(:,:,jp_sal), ldxios = lrxios ) ! before salt content isf_tsc trend 211 CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b' , risf_tsc_b(:,:,jp_tem), ldxios = lrxios ) ! before salt content isf_tsc trend 212 ELSE 213 fwfisf_b(:,:) = fwfisf(:,:) 214 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 215 ENDIF 216 ENDIF 59 IF ( ln_isfcav_mlt ) THEN 60 ! 61 ! before time step 62 IF ( kt /= nit000 ) THEN 63 risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:) 64 fwfisf_cav_b(:,:) = fwfisf_cav(:,:) 65 END IF 66 ! 67 ! compute tbl lvl/h 68 CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav) 69 ! 70 ! compute ice shelf melt 71 CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav) 72 ! 73 ! write restart variables (qoceisf, qhcisf, fwfisf for now and before) 74 IF (lrst_oce) CALL isfrst_write(kt, 'cav', risf_cav_tsc, fwfisf_cav) 75 ! 76 END IF 217 77 ! 218 IF( lrst_oce ) THEN 219 IF(lwp) WRITE(numout,*) 220 IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ', & 221 & 'at it= ', kt,' date= ', ndastp 222 IF(lwp) WRITE(numout,*) '~~~~' 223 IF( lwxios ) CALL iom_swap( cwxios_context ) 224 CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) , ldxios = lwxios ) 225 CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem), ldxios = lwxios ) 226 CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal), ldxios = lwxios ) 227 IF( lwxios ) CALL iom_swap( cxios_context ) 228 ENDIF 229 ! 230 END SUBROUTINE sbc_isf 231 232 233 INTEGER FUNCTION sbc_isf_alloc() 234 !!---------------------------------------------------------------------- 235 !! *** FUNCTION sbc_isf_rnf_alloc *** 236 !!---------------------------------------------------------------------- 237 sbc_isf_alloc = 0 ! set to zero if no array to be allocated 238 IF( .NOT. ALLOCATED( qisf ) ) THEN 239 ALLOCATE( risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj) , & 240 & rhisf_tbl(jpi,jpj) , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj) , & 241 & ttbl(jpi,jpj) , stbl(jpi,jpj) , utbl(jpi,jpj) , & 242 & vtbl(jpi, jpj) , risfLeff(jpi,jpj) , rhisf_tbl_0(jpi,jpj), & 243 & ralpha(jpi,jpj) , misfkt(jpi,jpj) , misfkb(jpi,jpj) , & 244 & STAT= sbc_isf_alloc ) 245 ! 246 CALL mpp_sum ( 'sbcisf', sbc_isf_alloc ) 247 IF( sbc_isf_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf_alloc: failed to allocate arrays.' ) 248 ! 249 ENDIF 250 END FUNCTION 251 252 253 SUBROUTINE sbc_isf_init 78 IF ( ln_isfpar_mlt ) THEN 79 ! 80 ! before time step 81 IF ( kt /= nit000 ) THEN 82 risf_par_tsc_b (:,:) = risf_par_tsc (:,:) 83 fwfisf_par_b(:,:) = fwfisf_par(:,:) 84 END IF 85 ! 86 ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl) 87 CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par) 88 ! 89 ! compute ice shelf melt 90 CALL isf_par( kt, risf_par_tsc, fwfisf_par) 91 ! 92 ! write restart variables (qoceisf, qhcisf, fwfisf for now and before) 93 IF (lrst_oce) CALL isfrst_write(kt, 'par', risf_par_tsc, fwfisf_par) 94 ! 95 END IF 96 ! 97 END SUBROUTINE isf_mlt 98 99 SUBROUTINE isf_mlt_init 254 100 !!--------------------------------------------------------------------- 255 !! *** ROUTINE sbc_isf_init *** 256 !! 257 !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 258 !! 259 !! ** Method : 4 parameterizations are available according to nn_isf 260 !! nn_isf = 1 : Realistic ice_shelf formulation 261 !! 2 : Beckmann & Goose parameterization 262 !! 3 : Specified runoff in deptht (Mathiot & al. ) 263 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 264 !!---------------------------------------------------------------------- 265 INTEGER :: ji, jj, jk ! loop index 266 INTEGER :: ik ! current level index 267 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 101 !! *** ROUTINE isfmlt_init *** 102 !! 103 !! ** Purpose : 104 !! 105 !! ** Method : 106 !! 107 !!---------------------------------------------------------------------- 268 108 INTEGER :: inum, ierror 269 109 INTEGER :: ios ! Local integer output status for namelist read 270 REAL(wp) :: zhk 271 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file 272 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 273 !!---------------------------------------------------------------------- 274 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 275 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 110 INTEGER :: ikt, ikb 111 INTEGER :: ji, jj 112 !!---------------------------------------------------------------------- 113 NAMELIST/namisf/ ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, sn_isfcav_fwf, & 114 & ln_isfpar_mlt, cn_isfpar_mlt, sn_isfpar_fwf, sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff 276 115 !!---------------------------------------------------------------------- 277 116 278 117 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 279 READ ( numnam_ref, nam sbc_isf, IOSTAT = ios, ERR = 901)280 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam sbc_isf in reference namelist', lwp )118 READ ( numnam_ref, namisf, IOSTAT = ios, ERR = 901) 119 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namisf in reference namelist', lwp ) 281 120 282 121 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 283 READ ( numnam_cfg, nam sbc_isf, IOSTAT = ios, ERR = 902 )284 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam sbc_isf in configuration namelist', lwp )285 IF(lwm) WRITE ( numond, nam sbc_isf )286 122 READ ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 ) 123 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namisf in configuration namelist', lwp ) 124 IF(lwm) WRITE ( numond, namisf ) 125 ! 287 126 IF(lwp) WRITE(numout,*) 288 IF(lwp) WRITE(numout,*) ' sbc_isf_init : heat flux of the ice shelf'127 IF(lwp) WRITE(numout,*) 'isf_init : ice shelf initialisation' 289 128 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 290 IF(lwp) WRITE(numout,*) ' Namelist namsbc_isf :' 291 IF(lwp) WRITE(numout,*) ' type ice shelf melting/freezing nn_isf = ', nn_isf 292 IF(lwp) WRITE(numout,*) ' bulk formulation (nn_isf=1 only) nn_isfblk = ', nn_isfblk 293 IF(lwp) WRITE(numout,*) ' thickness of the top boundary layer rn_hisf_tbl = ', rn_hisf_tbl 294 IF(lwp) WRITE(numout,*) ' gamma formulation nn_gammablk = ', nn_gammablk 295 IF(lwp) WRITE(numout,*) ' gammat coefficient rn_gammat0 = ', rn_gammat0 296 IF(lwp) WRITE(numout,*) ' gammas coefficient rn_gammas0 = ', rn_gammas0 297 IF(lwp) WRITE(numout,*) ' top drag coef. used (from namdrg_top) rn_Cd0 = ', r_Cdmin_top 298 299 300 ! 1 = presence of ISF 2 = bg03 parametrisation 301 ! 3 = rnf file for isf 4 = ISF fwf specified 302 ! option 1 and 4 need ln_isfcav = .true. (domzgr) 303 ! 304 ! Allocate public variable 305 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 129 IF(lwp) WRITE(numout,*) ' Namelist namisf :' 130 IF(lwp) WRITE(numout,*) ' melt inside the cavity ln_isfcav_mlt = ', ln_isfcav_mlt 131 IF ( ln_isfcav ) THEN 132 IF(lwp) WRITE(numout,*) ' melt formulation cn_isfcav_mlt = ', cn_isfcav_mlt 133 IF(lwp) WRITE(numout,*) ' thickness of the top boundary layer rn_htbl = ', rn_htbl 134 IF(lwp) WRITE(numout,*) ' gamma formulation cn_gammablk = ', cn_gammablk 135 IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN 136 IF(lwp) WRITE(numout,*) ' gammat coefficient rn_gammat0 = ', rn_gammat0 137 IF(lwp) WRITE(numout,*) ' gammas coefficient rn_gammas0 = ', rn_gammas0 138 IF(lwp) WRITE(numout,*) ' top drag coef. used (from namdrg_top) rn_Cd0 = ', r_Cdmin_top 139 END IF 140 END IF 141 IF(lwp) WRITE(numout,*) '' 142 IF(lwp) WRITE(numout,*) ' ice shelf melt parametrisation ln_isfpar_mlt = ', ln_isfpar_mlt 143 IF ( ln_isfpar_mlt ) THEN 144 IF(lwp) WRITE(numout,*) ' isf parametrisation formulation cn_isfpar_mlt = ', cn_isfpar_mlt 145 END IF 146 IF(lwp) WRITE(numout,*) '' 147 ! 148 ! Allocate public array 149 CALL isf_alloc() 306 150 ! 307 151 ! initialisation 308 qisf (:,:) = 0._wp ; fwfisf (:,:) = 0._wp 309 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp 310 ! 311 ! define isf tbl tickness, top and bottom indice 312 SELECT CASE ( nn_isf ) 313 CASE ( 1 ) 314 IF(lwp) WRITE(numout,*) 315 IF(lwp) WRITE(numout,*) ' ==>>> presence of under iceshelf seas (nn_isf = 1)' 316 rhisf_tbl(:,:) = rn_hisf_tbl 317 misfkt (:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 318 ! 319 CASE ( 2 , 3 ) 320 IF( .NOT.l_isfcpl ) THEN 321 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 322 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 323 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 324 ENDIF 325 ! read effective lenght (BG03) 326 IF( nn_isf == 2 ) THEN 327 IF(lwp) WRITE(numout,*) 328 IF(lwp) WRITE(numout,*) ' ==>>> bg03 parametrisation (nn_isf = 2)' 329 CALL iom_open( sn_Leff_isf%clname, inum ) 330 cvarLeff = TRIM(sn_Leff_isf%clvar) 331 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 332 CALL iom_close(inum) 333 ! 334 risfLeff = risfLeff*1000.0_wp !: convertion in m 335 ELSE 336 IF(lwp) WRITE(numout,*) 337 IF(lwp) WRITE(numout,*) ' ==>>> rnf file for isf (nn_isf = 3)' 338 ENDIF 339 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 340 CALL iom_open( sn_depmax_isf%clname, inum ) 341 cvarhisf = TRIM(sn_depmax_isf%clvar) 342 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 343 CALL iom_close(inum) 344 ! 345 CALL iom_open( sn_depmin_isf%clname, inum ) 346 cvarzisf = TRIM(sn_depmin_isf%clvar) 347 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 348 CALL iom_close(inum) 349 ! 350 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer 351 352 !! compute first level of the top boundary layer 353 DO ji = 1, jpi 354 DO jj = 1, jpj 355 ik = 2 356 !!gm potential bug: use gdepw_0 not _n 357 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO 358 misfkt(ji,jj) = ik-1 359 END DO 360 END DO 361 ! 362 CASE ( 4 ) 363 IF(lwp) WRITE(numout,*) 364 IF(lwp) WRITE(numout,*) ' ==>>> specified fresh water flux in ISF (nn_isf = 4)' 365 ! as in nn_isf == 1 366 rhisf_tbl(:,:) = rn_hisf_tbl 367 misfkt (:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 368 ! 369 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 370 IF( .NOT.l_isfcpl ) THEN 371 ALLOCATE( sf_fwfisf(1), STAT=ierror ) 372 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 373 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 374 ENDIF 375 ! 376 CASE DEFAULT 377 CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' ) 378 END SELECT 379 380 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 381 382 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 152 ! cav 153 misfkt_cav(:,:) = mikt(:,:) ; misfkb_cav(:,:) = mikt(:,:) 154 fwfisf_cav(:,:) = 0.0_wp ; fwfisf_cav_b(:,:) = 0.0_wp 155 rhisf_tbl_cav(:,:) = 1e-20 ; rfrac_tbl_cav(:,:) = 0.0_wp 156 risf_cav_tsc(:,:,:) = 0.0_wp ; risf_cav_tsc_b(:,:,:) = 0.0_wp 157 ! 158 mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:) 159 ! 160 ! par 161 misfkt_par(:,:) = 1 ; misfkb_par(:,:) = 1 162 fwfisf_par(:,:) = 0.0_wp ; fwfisf_par_b(:,:) = 0.0_wp 163 rhisf_tbl_par(:,:) = 1e-20 ; rfrac_tbl_par(:,:) = 0.0_wp 164 risf_par_tsc(:,:,:) = 0.0_wp ; risf_par_tsc_b(:,:,:) = 0.0_wp 165 ! 166 mskisf_par(:,:) = 0 167 ! 168 ! initialisation useful variable 169 r1_Lfusisf = 1._wp / rLfusisf 170 ! 383 171 DO jj = 1,jpj 384 172 DO ji = 1,jpi 385 ikt = misfkt(ji,jj) 386 ikb = misfkt(ji,jj) 387 ! thickness of boundary layer at least the top level thickness 388 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 389 390 ! determine the deepest level influenced by the boundary layer 391 DO jk = ikt+1, mbkt(ji,jj) 392 IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 393 END DO 394 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 395 misfkb(ji,jj) = ikb ! last wet level of the tbl 396 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 397 398 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 399 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 173 ikt = mikt(ji,jj) 174 ikb = mbkt(ji,jj) 175 bathy (ji,jj) = gdepw_0(ji,jj,ikb+1) 176 risfdep(ji,jj) = gdepw_0(ji,jj,ikt ) 400 177 END DO 401 178 END DO 402 403 IF( lwxios ) THEN 404 CALL iom_set_rstw_var_active('fwf_isf_b') 405 CALL iom_set_rstw_var_active('isf_hc_b') 406 CALL iom_set_rstw_var_active('isf_sc_b') 407 ENDIF 408 409 410 END SUBROUTINE sbc_isf_init 411 412 413 SUBROUTINE sbc_isf_bg03(kt) 414 !!--------------------------------------------------------------------- 415 !! *** ROUTINE sbc_isf_bg03 *** 416 !! 417 !! ** Purpose : add net heat and fresh water flux from ice shelf melting 418 !! into the adjacent ocean 419 !! 420 !! ** Method : See reference 421 !! 422 !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 423 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 424 !! (hereafter BG) 425 !! History : 06-02 (C. Wang) Original code 426 !!---------------------------------------------------------------------- 427 INTEGER, INTENT ( in ) :: kt 428 ! 429 INTEGER :: ji, jj, jk ! dummy loop index 430 INTEGER :: ik ! current level 431 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 432 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 433 REAL(wp) :: zt_frz ! freezing point temperature at depth z 434 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 435 !!---------------------------------------------------------------------- 436 ! 437 DO ji = 1, jpi 438 DO jj = 1, jpj 439 ik = misfkt(ji,jj) 440 !! Initialize arrays to 0 (each step) 441 zt_sum = 0.e0_wp 442 IF ( ik > 1 ) THEN 443 ! 1. -----------the average temperature between 200m and 600m --------------------- 444 DO jk = misfkt(ji,jj),misfkb(ji,jj) 445 ! Calculate freezing temperature 446 zpress = grav*rau0*gdept_n(ji,jj,ik)*1.e-04 447 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 448 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 449 END DO 450 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 451 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 452 ! For those corresponding to zonal boundary 453 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 454 & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 455 456 fwfisf(ji,jj) = qisf(ji,jj) / rLfusisf !fresh water flux kg/(m2s) 457 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 458 !add to salinity trend 459 ELSE 460 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 461 END IF 462 END DO 463 END DO 464 ! 465 END SUBROUTINE sbc_isf_bg03 466 467 468 SUBROUTINE sbc_isf_cav( kt ) 469 !!--------------------------------------------------------------------- 470 !! *** ROUTINE sbc_isf_cav *** 471 !! 472 !! ** Purpose : handle surface boundary condition under ice shelf 473 !! 474 !! ** Method : - 475 !! 476 !! ** Action : utau, vtau : remain unchanged 477 !! taum, wndm : remain unchanged 478 !! qns : update heat flux below ice shelf 479 !! emp, emps : update freshwater flux below ice shelf 480 !!--------------------------------------------------------------------- 481 INTEGER, INTENT(in) :: kt ! ocean time step 482 ! 483 INTEGER :: ji, jj ! dummy loop indices 484 INTEGER :: nit 485 LOGICAL :: lit 486 REAL(wp) :: zlamb1, zlamb2, zlamb3 487 REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 488 REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac 489 REAL(wp) :: zeps = 1.e-20_wp 490 REAL(wp) :: zerr 491 REAL(wp), DIMENSION(jpi,jpj) :: zfrz 492 REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas 493 REAL(wp), DIMENSION(jpi,jpj) :: zfwflx, zhtflx, zhtflx_b 494 !!--------------------------------------------------------------------- 495 ! 496 ! coeficient for linearisation of potential tfreez 497 ! Crude approximation for pressure (but commonly used) 498 IF ( l_useCT ) THEN ! linearisation from Jourdain et al. (2017) 499 zlamb1 =-0.0564_wp 500 zlamb2 = 0.0773_wp 501 zlamb3 =-7.8633e-8 * grav * rau0 502 ELSE ! linearisation from table 4 (Asay-Davis et al., 2015) 503 zlamb1 =-0.0573_wp 504 zlamb2 = 0.0832_wp 505 zlamb3 =-7.53e-8 * grav * rau0 506 ENDIF 507 ! 508 ! initialisation 509 zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 510 zhtflx (:,:) = 0.0_wp ; zhtflx_b(:,:) = 0.0_wp 511 zfwflx (:,:) = 0.0_wp 512 513 ! compute ice shelf melting 514 nit = 1 ; lit = .TRUE. 515 DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine 516 SELECT CASE ( nn_isfblk ) 517 CASE ( 1 ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 518 ! Calculate freezing temperature 519 CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 520 521 ! compute gammat every where (2d) 522 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 523 524 ! compute upward heat flux zhtflx and upward water flux zwflx 525 DO jj = 1, jpj 526 DO ji = 1, jpi 527 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 528 zfwflx(ji,jj) = - zhtflx(ji,jj)/rLfusisf 529 END DO 530 END DO 531 532 ! Compute heat flux and upward fresh water flux 533 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 534 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 535 536 CASE ( 2 ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 537 ! compute gammat every where (2d) 538 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 539 540 ! compute upward heat flux zhtflx and upward water flux zwflx 541 ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 542 DO jj = 1, jpj 543 DO ji = 1, jpi 544 ! compute coeficient to solve the 2nd order equation 545 zeps1 = rcp*rau0*zgammat(ji,jj) 546 zeps2 = rLfusisf*rau0*zgammas(ji,jj) 547 zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps) 548 zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 549 zeps6 = zeps4-ttbl(ji,jj) 550 zeps7 = zeps4-tsurf 551 zaqe = zlamb1 * (zeps1 + zeps3) 552 zaqer = 0.5_wp/MIN(zaqe,-zeps) 553 zbqe = zeps1*zeps6+zeps3*zeps7-zeps2 554 zcqe = zeps2*stbl(ji,jj) 555 zdis = zbqe*zbqe-4.0_wp*zaqe*zcqe 556 557 ! Presumably zdis can never be negative because gammas is very small compared to gammat 558 ! compute s freeze 559 zsfrz=(-zbqe-SQRT(zdis))*zaqer 560 IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 561 562 ! compute t freeze (eq. 22) 563 zfrz(ji,jj)=zeps4+zlamb1*zsfrz 564 565 ! zfwflx is upward water flux 566 ! zhtflx is upward heat flux (out of ocean) 567 ! compute the upward water and heat flux (eq. 28 and eq. 29) 568 zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 569 zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 570 END DO 571 END DO 572 573 ! compute heat and water flux 574 qisf (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 575 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 576 577 END SELECT 578 579 ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 580 IF ( nn_gammablk < 2 ) THEN ; lit = .FALSE. 581 ELSE 582 ! check total number of iteration 583 IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 584 ELSE ; nit = nit + 1 585 END IF 586 587 ! compute error between 2 iterations 588 ! if needed save gammat and compute zhtflx_b for next iteration 589 zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 590 IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 591 ELSE ; zhtflx_b(:,:) = zhtflx(:,:) 592 END IF 593 END IF 594 END DO 595 ! 596 CALL iom_put('isfgammat', zgammat) 597 CALL iom_put('isfgammas', zgammas) 598 ! 599 END SUBROUTINE sbc_isf_cav 600 601 602 SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 603 !!---------------------------------------------------------------------- 604 !! ** Purpose : compute the coefficient echange for heat flux 605 !! 606 !! ** Method : gamma assume constant or depends of u* and stability 607 !! 608 !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 609 !! Jenkins et al., 2010, JPO, p2298-2312 610 !!--------------------------------------------------------------------- 611 REAL(wp), DIMENSION(:,:), INTENT( out) :: pgt , pgs ! 612 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf ! 613 ! 614 INTEGER :: ji, jj ! loop index 615 INTEGER :: ikt ! local integer 616 REAL(wp) :: zdku, zdkv ! U, V shear 617 REAL(wp) :: zPr, zSc, zRc ! Prandtl, Scmidth and Richardson number 618 REAL(wp) :: zmob, zmols ! Monin Obukov length, coriolis factor at T point 619 REAL(wp) :: zbuofdep, zhnu ! Bouyancy length scale, sublayer tickness 620 REAL(wp) :: zhmax ! limitation of mol 621 REAL(wp) :: zetastar ! stability parameter 622 REAL(wp) :: zgmolet, zgmoles, zgturb ! contribution of modelecular sublayer and turbulence 623 REAL(wp) :: zcoef ! temporary coef 624 REAL(wp) :: zdep 625 REAL(wp) :: zeps = 1.0e-20_wp 626 REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant 627 REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 628 REAL(wp), DIMENSION(2) :: zts, zab 629 REAL(wp), DIMENSION(jpi,jpj) :: zustar ! U, V at T point and friction velocity 630 !!--------------------------------------------------------------------- 631 ! 632 SELECT CASE ( nn_gammablk ) 633 CASE ( 0 ) ! gamma is constant (specified in namelist) 634 !! ISOMIP formulation (Hunter et al, 2006) 635 pgt(:,:) = rn_gammat0 636 pgs(:,:) = rn_gammas0 637 638 CASE ( 1 ) ! gamma is assume to be proportional to u* 639 !! Jenkins et al., 2010, JPO, p2298-2312 640 !! Adopted by Asay-Davis et al. (2015) 641 !! compute ustar (eq. 24) 642 !!gm NB use pCdU here so that it will incorporate local boost of Cd0 and log layer case : 643 !! zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 644 !! or better : compute ustar in zdfdrg and use it here as well as in TKE, GLS and Co 645 !! 646 !! ===>>>> GM to be done this chrismas 647 !! 648 !!gm end 649 zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 650 651 !! Compute gammats 652 pgt(:,:) = zustar(:,:) * rn_gammat0 653 pgs(:,:) = zustar(:,:) * rn_gammas0 654 655 CASE ( 2 ) ! gamma depends of stability of boundary layer 656 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 657 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 658 !! compute ustar 659 zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 660 661 !! compute Pr and Sc number (can be improved) 662 zPr = 13.8_wp 663 zSc = 2432.0_wp 664 665 !! compute gamma mole 666 zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 667 zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 668 669 !! compute gamma 670 DO ji = 2, jpi 671 DO jj = 2, jpj 672 ikt = mikt(ji,jj) 673 674 IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think 675 pgt = rn_gammat0 676 pgs = rn_gammas0 677 ELSE 678 !! compute Rc number (as done in zdfric.F90) 679 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation 680 !!gm moreover, use Max(rn2,0) to take care of static instabilities.... 681 zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1) 682 ! ! shear of horizontal velocity 683 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & 684 & -un(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) ) 685 zdkv = zcoef * ( vn(ji ,jj-1,ikt ) + vn(ji,jj,ikt ) & 686 & -vn(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) ) 687 ! ! richardson number (minimum value set to zero) 688 zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 689 690 !! compute bouyancy 691 zts(jp_tem) = ttbl(ji,jj) 692 zts(jp_sal) = stbl(ji,jj) 693 zdep = gdepw_n(ji,jj,ikt) 694 ! 695 CALL eos_rab( zts, zdep, zab ) 696 ! 697 !! compute length scale 698 zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 699 700 !! compute Monin Obukov Length 701 ! Maximum boundary layer depth 702 zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp 703 ! Compute Monin obukhov length scale at the surface and Ekman depth: 704 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 705 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 706 707 !! compute eta* (stability parameter) 708 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp))) 709 710 !! compute the sublayer thickness 711 zhnu = 5 * znu / zustar(ji,jj) 712 713 !! compute gamma turb 714 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) & 715 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 716 717 !! compute gammats 718 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 719 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 720 END IF 721 END DO 722 END DO 723 CALL lbc_lnk_multi( 'sbcisf', pgt, 'T', 1., pgs, 'T', 1.) 724 END SELECT 725 ! 726 END SUBROUTINE sbc_isf_gammats 727 728 729 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 730 !!---------------------------------------------------------------------- 731 !! *** SUBROUTINE sbc_isf_tbl *** 732 !! 733 !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 734 !! 735 !!---------------------------------------------------------------------- 736 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pvarin 737 REAL(wp), DIMENSION(:,:) , INTENT( out) :: pvarout 738 CHARACTER(len=1), INTENT(in ) :: cd_ptin ! point of variable in/out 739 ! 740 INTEGER :: ji, jj, jk ! loop index 741 INTEGER :: ikt, ikb ! top and bottom index of the tbl 742 REAL(wp) :: ze3, zhk 743 REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl 744 !!---------------------------------------------------------------------- 745 746 ! initialisation 747 pvarout(:,:)=0._wp 748 749 SELECT CASE ( cd_ptin ) 750 CASE ( 'U' ) ! compute U in the top boundary layer at T- point 751 DO jj = 1,jpj 752 DO ji = 1,jpi 753 ikt = miku(ji,jj) ; ikb = miku(ji,jj) 754 ! thickness of boundary layer at least the top level thickness 755 zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) ) 756 757 ! determine the deepest level influenced by the boundary layer 758 DO jk = ikt+1, mbku(ji,jj) 759 IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 760 END DO 761 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 762 763 ! level fully include in the ice shelf boundary layer 764 DO jk = ikt, ikb - 1 765 ze3 = e3u_n(ji,jj,jk) 766 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 767 END DO 768 769 ! level partially include in ice shelf boundary layer 770 zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 771 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 772 END DO 773 END DO 774 DO jj = 2, jpj 775 DO ji = 2, jpi 776 !!gm a wet-point only average should be used here !!! 777 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 778 END DO 779 END DO 780 CALL lbc_lnk('sbcisf', pvarout,'T',-1.) 781 782 CASE ( 'V' ) ! compute V in the top boundary layer at T- point 783 DO jj = 1,jpj 784 DO ji = 1,jpi 785 ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 786 ! thickness of boundary layer at least the top level thickness 787 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt)) 788 789 ! determine the deepest level influenced by the boundary layer 790 DO jk = ikt+1, mbkv(ji,jj) 791 IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 792 END DO 793 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 794 795 ! level fully include in the ice shelf boundary layer 796 DO jk = ikt, ikb - 1 797 ze3 = e3v_n(ji,jj,jk) 798 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 799 END DO 800 801 ! level partially include in ice shelf boundary layer 802 zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 803 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 804 END DO 805 END DO 806 DO jj = 2, jpj 807 DO ji = 2, jpi 808 !!gm a wet-point only average should be used here !!! 809 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 810 END DO 811 END DO 812 CALL lbc_lnk('sbcisf', pvarout,'T',-1.) 813 814 CASE ( 'T' ) ! compute T in the top boundary layer at T- point 815 DO jj = 1,jpj 816 DO ji = 1,jpi 817 ikt = misfkt(ji,jj) 818 ikb = misfkb(ji,jj) 819 820 ! level fully include in the ice shelf boundary layer 821 DO jk = ikt, ikb - 1 822 ze3 = e3t_n(ji,jj,jk) 823 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 824 END DO 825 826 ! level partially include in ice shelf boundary layer 827 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 828 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 829 END DO 830 END DO 831 END SELECT 832 ! 833 ! mask mean tbl value 834 pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 835 ! 836 END SUBROUTINE sbc_isf_tbl 837 838 839 SUBROUTINE sbc_isf_div( phdivn ) 840 !!---------------------------------------------------------------------- 841 !! *** SUBROUTINE sbc_isf_div *** 842 !! 843 !! ** Purpose : update the horizontal divergence with the runoff inflow 844 !! 845 !! ** Method : 846 !! CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the 847 !! divergence and expressed in m/s 848 !! 849 !! ** Action : phdivn decreased by the runoff inflow 850 !!---------------------------------------------------------------------- 851 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdivn ! horizontal divergence 852 ! 853 INTEGER :: ji, jj, jk ! dummy loop indices 854 INTEGER :: ikt, ikb 855 REAL(wp) :: zhk 856 REAL(wp) :: zfact ! local scalar 857 !!---------------------------------------------------------------------- 858 ! 859 zfact = 0.5_wp 860 ! 861 IF(.NOT.ln_linssh ) THEN ! need to re compute level distribution of isf fresh water 862 DO jj = 1,jpj 863 DO ji = 1,jpi 864 ikt = misfkt(ji,jj) 865 ikb = misfkt(ji,jj) 866 ! thickness of boundary layer at least the top level thickness 867 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 868 869 ! determine the deepest level influenced by the boundary layer 870 DO jk = ikt, mbkt(ji,jj) 871 IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 872 END DO 873 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 874 misfkb(ji,jj) = ikb ! last wet level of the tbl 875 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 876 877 zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 878 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 879 END DO 880 END DO 881 END IF 882 ! 883 !== ice shelf melting distributed over several levels ==! 884 DO jj = 1,jpj 885 DO ji = 1,jpi 886 ikt = misfkt(ji,jj) 887 ikb = misfkb(ji,jj) 888 ! level fully include in the ice shelf boundary layer 889 DO jk = ikt, ikb - 1 890 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 891 & * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 892 END DO 893 ! level partially include in ice shelf boundary layer 894 phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 895 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 896 END DO 897 END DO 898 ! 899 END SUBROUTINE sbc_isf_div 900 179 ! 180 IF ( ln_isfcav_mlt ) THEN 181 ! 182 ! initialisation of cav variable 183 CALL isf_cav_init() 184 ! 185 ! read cav variable from restart 186 IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b) 187 ! 188 END IF 189 ! 190 IF ( ln_isfpar_mlt ) THEN 191 ! 192 ! initialisation of par variable 193 CALL isf_par_init() 194 ! 195 ! read par variable from restart 196 IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b) 197 ! 198 END IF 199 ! 200 END SUBROUTINE isf_mlt_init 901 201 !!====================================================================== 902 END MODULE sbcisf202 END MODULE isfmlt
Note: See TracChangeset
for help on using the changeset viewer.
