- Timestamp:
- 2015-11-29T20:28:41+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r5905 r5944 56 56 !: (first wet level and last level include in the tbl) 57 57 #else 58 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base58 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 59 59 #endif 60 60 61 61 62 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ?63 REAL(wp), PUBLIC, SAVE :: kappa= 1.54e-6_wp ! phycst ?64 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! phycst ?65 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! phycst ?66 REAL(wp), PUBLIC, SAVE :: lfusisf= 0.334e6_wp ! phycst ?62 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! phycst ? 63 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp ! phycst ? 64 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! phycst ? 65 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! phycst ? 66 REAL(wp), PUBLIC, SAVE :: rlfusisf = 0.334e6_wp ! phycst ? 67 67 68 68 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 69 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 70 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf file to be read 71 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 72 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the runoff file to be read 73 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 74 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf, sn_depmin_isf, sn_Leff_isf !: information about the runoff file to be read 69 CHARACTER(len=100), PUBLIC :: cn_dirisf = './' !: Root directory for location of ssr files 70 TYPE(FLD_N) , PUBLIC :: sn_fwfisf !: information about the isf file to be read 71 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 72 TYPE(FLD_N) , PUBLIC :: sn_rnfisf !: information about the runoff file to be read 73 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf 74 TYPE(FLD_N) , PUBLIC :: sn_depmax_isf !: information about the runoff file to be read ?? 75 TYPE(FLD_N) , PUBLIC :: sn_depmin_isf !: information about the runoff file to be read ?? 76 TYPE(FLD_N) , PUBLIC :: sn_Leff_isf !: information about the runoff file to be read ?? 75 77 76 78 !! * Substitutions … … 85 87 86 88 SUBROUTINE sbc_isf(kt) 87 INTEGER, INTENT(in) :: kt ! ocean time step88 INTEGER :: ji, jj, jk ! loop index89 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer90 INTEGER :: inum, ierror91 INTEGER :: ios ! Local integer output status for namelist read92 REAL(wp) :: zhk93 REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)94 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file95 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale96 !97 89 !!--------------------------------------------------------------------- 98 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf , & 99 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 100 ! 101 ! allocation 102 CALL wrk_alloc( jpi,jpj, zt_frz, zdep ) 90 !! *** ROUTINE sbc_isf *** 91 !! 92 !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 93 !! melting and freezing 94 !! 95 !! ** Method : 4 parameterizations are available according to nn_isf 96 !! nn_isf = 1 : Realistic ice_shelf formulation 97 !! 2 : Beckmann & Goose parameterization 98 !! 3 : Specified runoff in deptht (Mathiot & al. ) 99 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 100 !!---------------------------------------------------------------------- 101 INTEGER, INTENT( in ) :: kt ! ocean time step 102 ! 103 INTEGER :: ji, jj ! loop index 104 REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 105 !!--------------------------------------------------------------------- 103 106 ! 104 107 ! ! ====================== ! 105 108 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 106 109 ! ! ====================== ! 107 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 108 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 109 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 110 111 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 112 READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 113 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 114 IF(lwm) WRITE ( numond, namsbc_isf ) 115 116 117 IF ( lwp ) WRITE(numout,*) 118 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 119 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 120 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 121 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 122 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 123 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk 124 IF ( lwp ) WRITE(numout,*) ' rn_hisf_tbl = ', rn_hisf_tbl 125 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 126 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 127 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 128 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 110 CALL sbc_isf_init 111 ! ! ---------------------------------------- ! 112 ELSE ! Swap of forcing fields ! 113 ! ! ---------------------------------------- ! 114 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000 115 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine 129 116 ! 130 ! Allocate public variable131 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )132 !133 ! initialisation134 qisf(:,:) = 0._wp ; fwfisf (:,:) = 0._wp135 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp136 !137 ! define isf tbl tickness, top and bottom indice138 IF (nn_isf == 1) THEN139 rhisf_tbl(:,:) = rn_hisf_tbl140 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv141 ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN142 ALLOCATE( sf_rnfisf(1), STAT=ierror )143 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )144 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )145 146 !: read effective lenght (BG03)147 IF (nn_isf == 2) THEN148 cvarLeff = 'soLeff'149 CALL iom_open( sn_Leff_isf%clname, inum )150 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)151 CALL iom_close(inum)152 !153 risfLeff = risfLeff*1000.0_wp !: convertion in m154 END IF155 156 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)157 CALL iom_open( sn_depmax_isf%clname, inum )158 cvarhisf = TRIM(sn_depmax_isf%clvar)159 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base160 CALL iom_close(inum)161 !162 CALL iom_open( sn_depmin_isf%clname, inum )163 cvarzisf = TRIM(sn_depmin_isf%clvar)164 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base165 CALL iom_close(inum)166 !167 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer168 169 !! compute first level of the top boundary layer170 DO ji = 1, jpi171 DO jj = 1, jpj172 jk = 2173 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ; jk = jk + 1 ; END DO174 misfkt(ji,jj) = jk-1175 END DO176 END DO177 178 ELSE IF ( nn_isf == 4 ) THEN179 ! as in nn_isf == 1180 rhisf_tbl(:,:) = rn_hisf_tbl181 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv182 183 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)184 ALLOCATE( sf_fwfisf(1), STAT=ierror )185 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )186 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )187 END IF188 189 rhisf_tbl_0(:,:) = rhisf_tbl(:,:)190 191 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf192 DO jj = 1,jpj193 DO ji = 1,jpi194 ikt = misfkt(ji,jj)195 ikb = misfkt(ji,jj)196 ! thickness of boundary layer at least the top level thickness197 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))198 199 ! determine the deepest level influenced by the boundary layer200 DO jk = ikt+1, mbkt(ji,jj)201 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk202 END DO203 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.204 misfkb(ji,jj) = ikb ! last wet level of the tbl205 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)206 207 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1208 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer209 END DO210 END DO211 212 117 END IF 213 118 214 ! ! ---------------------------------------- !215 IF( kt /= nit000 ) THEN ! Swap of forcing fields !216 ! ! ---------------------------------------- !217 fwfisf_b (:,: ) = fwfisf (:,: ) ! Swap the ocean forcing fields except at nit000218 risf_tsc_b(:,:,:) = risf_tsc(:,:,:) ! where before fields are set at the end of the routine219 !220 ENDIF221 222 119 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 223 224 225 ! compute salf and heat flux 226 IF (nn_isf == 1) THEN 227 ! realistic ice shelf formulation 120 ! allocation 121 CALL wrk_alloc( jpi,jpj, zt_frz, zdep ) 122 123 ! compute salt and heat flux 124 SELECT CASE ( nn_isf ) 125 CASE ( 1 ) ! realistic ice shelf formulation 228 126 ! compute T/S/U/V for the top boundary layer 229 127 CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') … … 239 137 CALL sbc_isf_cav (kt) 240 138 241 ELSE IF (nn_isf == 2) THEN 242 ! Beckmann and Goosse parametrisation 139 CASE ( 2 ) ! Beckmann and Goosse parametrisation 243 140 stbl(:,:) = soce 244 141 CALL sbc_isf_bg03(kt) 245 142 246 ELSE IF (nn_isf == 3) THEN 247 ! specified runoff in depth (Mathiot et al., XXXX in preparation) 143 CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation) 248 144 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 249 145 fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 250 qisf(:,:) = fwfisf(:,:) * lfusisf! heat flux146 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 251 147 stbl(:,:) = soce 252 148 253 ELSE IF (nn_isf == 4) THEN 254 ! specified fwf and heat flux forcing beneath the ice shelf 149 CASE ( 4 ) ! specified fwf and heat flux forcing beneath the ice shelf 255 150 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 256 151 fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1) ! fwf flux from the isf (fwfisf <0 mean melting) 257 qisf(:,:) = fwfisf(:,:) * lfusisf! heat flux152 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 258 153 stbl(:,:) = soce 259 END IF 154 155 END SELECT 260 156 261 157 ! compute tsc due to isf … … 276 172 CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 277 173 CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 278 CALL lbc_lnk(fwfisf(:,:) ,'T',1.)279 CALL lbc_lnk(qisf(:,:) ,'T',1.)174 CALL lbc_lnk(fwfisf(:,:) ,'T',1.) 175 CALL lbc_lnk(qisf(:,:) ,'T',1.) 280 176 281 177 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! … … 290 186 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 291 187 END IF 292 END IF188 END IF 293 189 ! 294 190 ! output 295 IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf) 296 IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf) 191 ! JMM : iom_use not necessary here, qisf and fwfisf are always computed. 192 ! If not required in iodef.xml, iom_put does not do anything 193 ! IF( iom_use('qisf' ) ) CALL iom_put('qisf' , qisf) 194 ! IF( iom_use('fwfisf') ) CALL iom_put('fwfisf', fwfisf) 195 CALL iom_put('qisf' , qisf) 196 CALL iom_put('fwfisf', fwfisf) 197 198 ! deallocation 199 CALL wrk_dealloc( jpi,jpj, zt_frz, zdep ) 297 200 END IF 298 299 ! deallocation300 CALL wrk_dealloc( jpi,jpj, zt_frz, zdep )301 201 302 202 END SUBROUTINE sbc_isf … … 315 215 & STAT= sbc_isf_alloc ) 316 216 ! 317 IF( lk_mpp 217 IF( lk_mpp ) CALL mpp_sum ( sbc_isf_alloc ) 318 218 IF( sbc_isf_alloc /= 0 ) CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 319 219 ! 320 END IF220 END IF 321 221 END FUNCTION 322 222 223 SUBROUTINE sbc_isf_init 224 !!--------------------------------------------------------------------- 225 !! *** ROUTINE sbc_isf_init *** 226 !! 227 !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 228 !! 229 !! ** Method : 4 parameterizations are available according to nn_isf 230 !! nn_isf = 1 : Realistic ice_shelf formulation 231 !! 2 : Beckmann & Goose parameterization 232 !! 3 : Specified runoff in deptht (Mathiot & al. ) 233 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 234 !!---------------------------------------------------------------------- 235 INTEGER :: ji, jj, jk ! loop index 236 INTEGER :: ik ! current level index 237 INTEGER :: ikt, ikb ! top and bottom level of the isf boundary layer 238 INTEGER :: inum, ierror 239 INTEGER :: ios ! Local integer output status for namelist read 240 REAL(wp) :: zhk 241 CHARACTER(len=256) :: cvarzisf, cvarhisf ! name for isf file 242 CHARACTER(LEN=32 ) :: cvarLeff ! variable name for efficient Length scale 243 !!---------------------------------------------------------------------- 244 NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 245 & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 246 !!---------------------------------------------------------------------- 247 248 REWIND( numnam_ref ) ! Namelist namsbc_rnf in reference namelist : Runoffs 249 READ ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 250 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 251 252 REWIND( numnam_cfg ) ! Namelist namsbc_rnf in configuration namelist : Runoffs 253 READ ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 254 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 255 IF(lwm) WRITE ( numond, namsbc_isf ) 256 257 IF ( lwp ) WRITE(numout,*) 258 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 259 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 260 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 261 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 262 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 263 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk 264 IF ( lwp ) WRITE(numout,*) ' rn_hisf_tbl = ', rn_hisf_tbl 265 IF ( lwp ) WRITE(numout,*) ' nn_gammablk = ', nn_gammablk 266 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 267 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 268 IF ( lwp ) WRITE(numout,*) ' rn_tfri2 = ', rn_tfri2 269 ! 270 ! Allocate public variable 271 IF ( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 272 ! 273 ! initialisation 274 qisf(:,:) = 0._wp ; fwfisf (:,:) = 0._wp 275 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp 276 ! 277 ! define isf tbl tickness, top and bottom indice 278 SELECT CASE ( nn_isf ) 279 CASE ( 1 ) 280 rhisf_tbl(:,:) = rn_hisf_tbl 281 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 282 283 CASE ( 2 , 3 ) 284 ALLOCATE( sf_rnfisf(1), STAT=ierror ) 285 ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 286 CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 287 288 ! read effective lenght (BG03) 289 IF (nn_isf == 2) THEN 290 CALL iom_open( sn_Leff_isf%clname, inum ) 291 cvarLeff = TRIM(sn_Leff_isf%clvar) 292 CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 293 CALL iom_close(inum) 294 ! 295 risfLeff = risfLeff*1000.0_wp !: convertion in m 296 END IF 297 298 ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 299 CALL iom_open( sn_depmax_isf%clname, inum ) 300 cvarhisf = TRIM(sn_depmax_isf%clvar) 301 CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 302 CALL iom_close(inum) 303 ! 304 CALL iom_open( sn_depmin_isf%clname, inum ) 305 cvarzisf = TRIM(sn_depmin_isf%clvar) 306 CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 307 CALL iom_close(inum) 308 ! 309 rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:) !: tickness isf boundary layer 310 311 !! compute first level of the top boundary layer 312 DO ji = 1, jpi 313 DO jj = 1, jpj 314 ik = 2 315 DO WHILE ( ik <= mbkt(ji,jj) .AND. fsdepw(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO 316 misfkt(ji,jj) = ik-1 317 END DO 318 END DO 319 320 CASE ( 4 ) 321 ! as in nn_isf == 1 322 rhisf_tbl(:,:) = rn_hisf_tbl 323 misfkt(:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 324 325 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 326 ALLOCATE( sf_fwfisf(1), STAT=ierror ) 327 ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 328 CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 329 330 END SELECT 331 332 rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 333 334 ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 335 DO jj = 1,jpj 336 DO ji = 1,jpi 337 ikt = misfkt(ji,jj) 338 ikb = misfkt(ji,jj) 339 ! thickness of boundary layer at least the top level thickness 340 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 341 342 ! determine the deepest level influenced by the boundary layer 343 DO jk = ikt+1, mbkt(ji,jj) 344 IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 345 END DO 346 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 347 misfkb(ji,jj) = ikb ! last wet level of the tbl 348 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 349 350 zhk = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 351 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer 352 END DO 353 END DO 354 355 END SUBROUTINE sbc_isf_init 356 323 357 SUBROUTINE sbc_isf_bg03(kt) 324 !!========================================================================== 325 !! *** SUBROUTINE sbcisf_bg03 *** 326 !! add net heat and fresh water flux from ice shelf melting 327 !! into the adjacent ocean using the parameterisation by 328 !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 329 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 330 !! (hereafter BG) 331 !!========================================================================== 332 !!---------------------------------------------------------------------- 333 !! sbc_isf_bg03 : routine called from sbcmod 334 !!---------------------------------------------------------------------- 335 !! 336 !! ** Purpose : Add heat and fresh water fluxes due to ice shelf melting 337 !! ** Reference : Beckmann et Goosse, 2003, Ocean Modelling 338 !! 339 !! History : 340 !! ! 06-02 (C. Wang) Original code 341 !!---------------------------------------------------------------------- 342 343 INTEGER, INTENT ( in ) :: kt 344 345 INTEGER :: ji, jj, jk !temporary integer 346 347 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 348 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 349 REAL(wp) :: zt_frz ! freezing point temperature at depth z 350 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 351 352 !!---------------------------------------------------------------------- 353 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 354 ! 355 356 ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run ) 357 DO ji = 1, jpi 358 DO jj = 1, jpj 359 jk = misfkt(ji,jj) 360 !! Initialize arrays to 0 (each step) 361 zt_sum = 0.e0_wp 362 IF ( jk .GT. 1 ) THEN 363 ! 1. -----------the average temperature between 200m and 600m --------------------- 364 DO jk = misfkt(ji,jj),misfkb(ji,jj) 365 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 366 ! after verif with UNESCO, wrong sign in BG eq. 2 367 ! Calculate freezing temperature 368 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 369 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 370 ENDDO 371 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 372 373 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 374 ! For those corresponding to zonal boundary 375 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 376 & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,jk) 358 !!--------------------------------------------------------------------- 359 !! *** ROUTINE sbc_isf_bg03 *** 360 !! 361 !! ** Purpose : add net heat and fresh water flux from ice shelf melting 362 !! into the adjacent ocean 363 !! 364 !! ** Method : See reference 365 !! 366 !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 367 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 368 !! (hereafter BG) 369 !! History : 370 !! 06-02 (C. Wang) Original code 371 !!---------------------------------------------------------------------- 372 INTEGER, INTENT ( in ) :: kt 373 ! 374 INTEGER :: ji, jj, jk ! dummy loop index 375 INTEGER :: ik ! current level 376 REAL(wp) :: zt_sum ! sum of the temperature between 200m and 600m 377 REAL(wp) :: zt_ave ! averaged temperature between 200m and 600m 378 REAL(wp) :: zt_frz ! freezing point temperature at depth z 379 REAL(wp) :: zpress ! pressure to compute the freezing point in depth 380 !!---------------------------------------------------------------------- 381 382 IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 383 ! 384 DO ji = 1, jpi 385 DO jj = 1, jpj 386 ik = misfkt(ji,jj) 387 !! Initialize arrays to 0 (each step) 388 zt_sum = 0.e0_wp 389 IF ( ik > 1 ) THEN 390 ! 1. -----------the average temperature between 200m and 600m --------------------- 391 DO jk = misfkt(ji,jj),misfkb(ji,jj) 392 ! freezing point temperature at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 393 ! after verif with UNESCO, wrong sign in BG eq. 2 394 ! Calculate freezing temperature 395 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 396 zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp 397 END DO 398 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 399 ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 400 ! For those corresponding to zonal boundary 401 qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave & 402 & * r1_e12t(ji,jj) * tmask(ji,jj,jk) 377 403 378 fwfisf(ji,jj) = qisf(ji,jj) / lfusisf !fresh water flux kg/(m2s) 379 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 380 !add to salinity trend 381 ELSE 382 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 383 END IF 384 ENDDO 385 ENDDO 386 ! 387 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03') 404 fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf !fresh water flux kg/(m2s) 405 fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 406 !add to salinity trend 407 ELSE 408 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 409 END IF 410 END DO 411 END DO 412 ! 413 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03') 414 388 415 END SUBROUTINE sbc_isf_bg03 389 416 390 417 SUBROUTINE sbc_isf_cav( kt ) 391 418 !!--------------------------------------------------------------------- 392 419 !! *** ROUTINE sbc_isf_cav *** … … 403 430 INTEGER, INTENT(in) :: kt ! ocean time step 404 431 ! 405 REAL(wp), DIMENSION(:,:), POINTER :: zfrz 406 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 407 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b 432 INTEGER :: ji, jj ! dummy loop indices 433 INTEGER :: nit 408 434 REAL(wp) :: zlamb1, zlamb2, zlamb3 409 435 REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 … … 411 437 REAL(wp) :: zeps = 1.e-20_wp 412 438 REAL(wp) :: zerr 413 INTEGER :: ji, jj ! dummy loop indices 414 INTEGER :: nit 439 REAL(wp), DIMENSION(:,:), POINTER :: zfrz 440 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 441 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b 415 442 LOGICAL :: lit 416 443 !!--------------------------------------------------------------------- 417 !418 444 ! coeficient for linearisation of potential tfreez 419 445 ! Crude approximation for pressure (but commonly used) 420 zlamb1 =-0.0573_wp421 zlamb2 = 0.0832_wp422 zlamb3 =-7.53e-08* grav * rau0446 zlamb1 =-0.0573_wp 447 zlamb2 = 0.0832_wp 448 zlamb3 =-7.53e-08_wp * grav * rau0 423 449 IF( nn_timing == 1 ) CALL timing_start('sbc_isf_cav') 424 450 ! … … 427 453 428 454 ! initialisation 429 zgammat(:,:) =rn_gammat0 ; zgammas (:,:)=rn_gammas0;430 zhtflx (:,:) =0.0_wp ; zhtflx_b(:,:)=0.0_wp ;431 zfwflx (:,:) =0.0_wp455 zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 456 zhtflx (:,:) = 0.0_wp ; zhtflx_b(:,:) = 0.0_wp 457 zfwflx (:,:) = 0.0_wp 432 458 433 459 ! compute ice shelf melting 434 460 nit = 1 ; lit = .TRUE. 435 461 DO WHILE ( lit ) ! maybe just a constant number of iteration as in blk_core is fine 436 IF (nn_isfblk == 1) THEN437 !ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)462 SELECT CASE ( nn_isfblk ) 463 CASE ( 1 ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 438 464 ! Calculate freezing temperature 439 465 CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) … … 446 472 DO ji = 1, jpi 447 473 zhtflx(ji,jj) = zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 448 zfwflx(ji,jj) = - zhtflx(ji,jj)/ lfusisf474 zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 449 475 END DO 450 476 END DO … … 454 480 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 455 481 456 ELSE IF (nn_isfblk == 2 ) THEN 457 ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 482 CASE ( 2 ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 458 483 ! compute gammat every where (2d) 459 484 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) … … 464 489 DO ji = 1, jpi 465 490 ! compute coeficient to solve the 2nd order equation 466 zeps1 =rcp*rau0*zgammat(ji,jj)467 zeps2 =lfusisf*rau0*zgammas(ji,jj)468 zeps3 =rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps)469 zeps4 =zlamb2+zlamb3*risfdep(ji,jj)470 zeps6 =zeps4-ttbl(ji,jj)471 zeps7 =zeps4-tsurf472 zaqe =zlamb1 * (zeps1 + zeps3)473 zaqer =0.5/MIN(zaqe,-zeps)474 zbqe =zeps1*zeps6+zeps3*zeps7-zeps2475 zcqe =zeps2*stbl(ji,jj)476 zdis =zbqe*zbqe-4.0*zaqe*zcqe491 zeps1 = rcp*rau0*zgammat(ji,jj) 492 zeps2 = rlfusisf*rau0*zgammas(ji,jj) 493 zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 494 zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 495 zeps6 = zeps4-ttbl(ji,jj) 496 zeps7 = zeps4-tsurf 497 zaqe = zlamb1 * (zeps1 + zeps3) 498 zaqer = 0.5_wp/MIN(zaqe,-zeps) 499 zbqe = zeps1*zeps6+zeps3*zeps7-zeps2 500 zcqe = zeps2*stbl(ji,jj) 501 zdis = zbqe*zbqe-4.0_wp*zaqe*zcqe 477 502 478 503 ! Presumably zdis can never be negative because gammas is very small compared to gammat 479 504 ! compute s freeze 480 505 zsfrz=(-zbqe-SQRT(zdis))*zaqer 481 IF ( zsfrz .LT.0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer506 IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 482 507 483 508 ! compute t freeze (eq. 22) … … 496 521 fwfisf(:,:) = zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 497 522 498 END IF523 END SELECT 499 524 500 525 ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 501 IF ( nn_gammablk .LT.2 ) THEN ; lit = .FALSE.526 IF ( nn_gammablk < 2 ) THEN ; lit = .FALSE. 502 527 ELSE 503 528 ! check total number of iteration 504 IF (nit .GE.100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )505 ELSE 506 END IF529 IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 530 ELSE ; nit = nit + 1 531 END IF 507 532 508 533 ! compute error between 2 iterations 509 534 ! if needed save gammat and compute zhtflx_b for next iteration 510 535 zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 511 IF ( zerr .LE. 0.01) THEN ; lit = .FALSE.512 ELSE ; zhtflx_b(:,:) = zhtflx(:,:)513 END IF536 IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 537 ELSE ; zhtflx_b(:,:) = zhtflx(:,:) 538 END IF 514 539 END IF 515 540 END DO 516 541 ! 517 ! output 518 IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat) 519 IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 542 ! output see JMM comment above 543 ! IF( iom_use('isfgammat') ) CALL iom_put('isfgammat', zgammat) 544 ! IF( iom_use('isfgammas') ) CALL iom_put('isfgammas', zgammas) 545 CALL iom_put('isfgammat', zgammat) 546 CALL iom_put('isfgammas', zgammas) 520 547 ! 521 548 CALL wrk_dealloc( jpi,jpj, zfrz , zgammat, zgammas ) … … 526 553 END SUBROUTINE sbc_isf_cav 527 554 528 SUBROUTINE sbc_isf_gammats( gt, gs, zqhisf, zqwisf )555 SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 529 556 !!---------------------------------------------------------------------- 530 557 !! ** Purpose : compute the coefficient echange for heat flux … … 535 562 !! Jenkins et al., 2010, JPO, p2298-2312 536 563 !!--------------------------------------------------------------------- 537 REAL(wp), DIMENSION(:,:), INTENT(out) :: gt,gs538 REAL(wp), DIMENSION(:,:), INTENT(in ) :: zqhisf, zqwisf539 564 REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs 565 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf 566 ! 540 567 INTEGER :: ikt 541 INTEGER :: ji, jj! loop index568 INTEGER :: ji, jj ! loop index 542 569 REAL(wp), DIMENSION(:,:), POINTER :: zustar ! U, V at T point and friction velocity 543 570 REAL(wp) :: zdku, zdkv ! U, V shear … … 551 578 REAL(wp) :: zdep 552 579 REAL(wp) :: zeps = 1.0e-20_wp 553 REAL(wp), PARAMETER :: zxsiN = 0.052 ! dimensionless constant554 REAL(wp), PARAMETER :: znu = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)580 REAL(wp), PARAMETER :: zxsiN = 0.052_wp ! dimensionless constant 581 REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 555 582 REAL(wp), DIMENSION(2) :: zts, zab 556 583 !!--------------------------------------------------------------------- 557 584 CALL wrk_alloc( jpi,jpj, zustar ) 558 585 ! 559 IF ( nn_gammablk == 0 ) THEN 560 !! gamma is constant (specified in namelist) 561 !! ISOMIP formulation (Hunter et al, 2006) 562 gt(:,:) = rn_gammat0 563 gs(:,:) = rn_gammas0 564 565 ELSE IF ( nn_gammablk == 1 ) THEN 566 !! gamma is assume to be proportional to u* 567 !! Jenkins et al., 2010, JPO, p2298-2312 568 !! Adopted by Asay-Davis et al. (2015) 569 570 !! compute ustar (eq. 24) 586 SELECT CASE ( nn_gammablk ) 587 CASE ( 0 ) ! gamma is constant (specified in namelist) 588 !! ISOMIP formulation (Hunter et al, 2006) 589 pgt(:,:) = rn_gammat0 590 pgs(:,:) = rn_gammas0 591 592 CASE ( 1 ) ! gamma is assume to be proportional to u* 593 !! Jenkins et al., 2010, JPO, p2298-2312 594 !! Adopted by Asay-Davis et al. (2015) 595 596 !! compute ustar (eq. 24) 571 597 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 572 598 573 !! Compute gammats574 gt(:,:) = zustar(:,:) * rn_gammat0575 gs(:,:) = zustar(:,:) * rn_gammas0599 !! Compute gammats 600 pgt(:,:) = zustar(:,:) * rn_gammat0 601 pgs(:,:) = zustar(:,:) * rn_gammas0 576 602 577 ELSE IF ( nn_gammablk == 2 ) THEN 578 !! gamma depends of stability of boundary layer 579 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 580 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 581 !! compute ustar 603 CASE ( 2 ) ! gamma depends of stability of boundary layer 604 !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 605 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 606 !! compute ustar 582 607 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 583 608 584 !! compute Pr and Sc number (can be improved)585 zPr = 13.8 586 zSc = 2432.0 587 588 !! compute gamma mole589 zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0590 zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0591 592 !! compute gamma609 !! compute Pr and Sc number (can be improved) 610 zPr = 13.8_wp 611 zSc = 2432.0_wp 612 613 !! compute gamma mole 614 zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 615 zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 616 617 !! compute gamma 593 618 DO ji=2,jpi 594 619 DO jj=2,jpj … … 596 621 597 622 IF (zustar(ji,jj) == 0._wp) THEN ! only for kt = 1 I think 598 gt = rn_gammat0599 gs = rn_gammas0623 pgt = rn_gammat0 624 pgs = rn_gammas0 600 625 ELSE 601 !! compute Rc number (as done in zdfric.F90)602 zcoef = 0.5 / fse3w(ji,jj,ikt)626 !! compute Rc number (as done in zdfric.F90) 627 zcoef = 0.5_wp / fse3w(ji,jj,ikt) 603 628 ! ! shear of horizontal velocity 604 629 zdku = zcoef * ( un(ji-1,jj ,ikt ) + un(ji,jj,ikt ) & … … 609 634 zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 610 635 611 !! compute bouyancy636 !! compute bouyancy 612 637 zts(jp_tem) = ttbl(ji,jj) 613 638 zts(jp_sal) = stbl(ji,jj) … … 616 641 CALL eos_rab( zts, zdep, zab ) 617 642 ! 618 !! compute length scale619 zbuofdep = grav * ( zab(jp_tem) * zqhisf(ji,jj) - zab(jp_sal) * zqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!620 621 !! compute Monin Obukov Length643 !! compute length scale 644 zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) ) !!!!!!!!!!!!!!!!!!!!!!!!!!!! 645 646 !! compute Monin Obukov Length 622 647 ! Maximum boundary layer depth 623 648 zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp … … 626 651 zmols = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 627 652 628 !! compute eta* (stability parameter)653 !! compute eta* (stability parameter) 629 654 zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 630 655 631 !! compute the sublayer thickness656 !! compute the sublayer thickness 632 657 zhnu = 5 * znu / zustar(ji,jj) 633 658 634 !! compute gamma turb659 !! compute gamma turb 635 660 zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 636 661 & + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 637 662 638 !! compute gammats639 gt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)640 gs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)663 !! compute gammats 664 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 665 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 641 666 END IF 642 667 END DO 643 668 END DO 644 CALL lbc_lnk( gt(:,:),'T',1.)645 CALL lbc_lnk( gs(:,:),'T',1.)646 END IF669 CALL lbc_lnk(pgt(:,:),'T',1.) 670 CALL lbc_lnk(pgs(:,:),'T',1.) 671 END SELECT 647 672 CALL wrk_dealloc( jpi,jpj, zustar ) 648 673 649 END SUBROUTINE 650 651 SUBROUTINE sbc_isf_tbl( varin, varout, cptin )674 END SUBROUTINE sbc_isf_gammats 675 676 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 652 677 !!---------------------------------------------------------------------- 653 678 !! *** SUBROUTINE sbc_isf_tbl *** … … 656 681 !! 657 682 !!---------------------------------------------------------------------- 658 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin 659 REAL(wp), DIMENSION(:,:) , INTENT(out):: varout 660 661 CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out 662 683 REAL(wp), DIMENSION(:,:,:), INTENT( in ) :: pvarin 684 REAL(wp), DIMENSION(:,:) , INTENT( out ) :: pvarout 685 CHARACTER(len=1), INTENT( in ) :: cd_ptin ! point of variable in/out 686 ! 663 687 REAL(wp) :: ze3, zhk 664 688 REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 665 689 666 INTEGER :: ji, jj,jk! loop index667 INTEGER :: ikt, ikb ! top and bottom index of the tbl668 690 INTEGER :: ji, jj, jk ! loop index 691 INTEGER :: ikt, ikb ! top and bottom index of the tbl 692 !!---------------------------------------------------------------------- 669 693 ! allocation 670 694 CALL wrk_alloc( jpi,jpj, zhisf_tbl) 671 695 672 696 ! initialisation 673 varout(:,:)=0._wp697 pvarout(:,:)=0._wp 674 698 675 ! compute U in the top boundary layer at T- point676 IF (cptin == 'U') THEN699 SELECT CASE ( cd_ptin ) 700 CASE ( 'U' ) ! compute U in the top boundary layer at T- point 677 701 DO jj = 1,jpj 678 702 DO ji = 1,jpi 679 703 ikt = miku(ji,jj) ; ikb = miku(ji,jj) 680 ! thickness of boundary layer at least the top level thickness704 ! thickness of boundary layer at least the top level thickness 681 705 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt)) 682 706 683 ! determine the deepest level influenced by the boundary layer707 ! determine the deepest level influenced by the boundary layer 684 708 DO jk = ikt+1, mbku(ji,jj) 685 IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) .LT.zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk709 IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 686 710 END DO 687 711 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 688 712 689 ! level fully include in the ice shelf boundary layer713 ! level fully include in the ice shelf boundary layer 690 714 DO jk = ikt, ikb - 1 691 715 ze3 = fse3u_n(ji,jj,jk) 692 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3693 END DO 694 695 ! level partially include in ice shelf boundary layer716 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 717 END DO 718 719 ! level partially include in ice shelf boundary layer 696 720 zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 697 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,ikb) * (1._wp - zhk)721 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 698 722 END DO 699 723 END DO 700 724 DO jj = 2,jpj 701 725 DO ji = 2,jpi 702 varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji-1,jj)) 703 END DO 704 END DO 705 CALL lbc_lnk(varout,'T',-1.) 706 END IF 707 708 ! compute V in the top boundary layer at T- point 709 IF (cptin == 'V') THEN 726 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 727 END DO 728 END DO 729 CALL lbc_lnk(pvarout,'T',-1.) 730 731 CASE ( 'V' ) ! compute V in the top boundary layer at T- point 710 732 DO jj = 1,jpj 711 733 DO ji = 1,jpi 712 734 ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 713 ! thickness of boundary layer at least the top level thickness735 ! thickness of boundary layer at least the top level thickness 714 736 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt)) 715 737 716 ! determine the deepest level influenced by the boundary layer738 ! determine the deepest level influenced by the boundary layer 717 739 DO jk = ikt+1, mbkv(ji,jj) 718 IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) .LT.zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk740 IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 719 741 END DO 720 742 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 721 743 722 ! level fully include in the ice shelf boundary layer744 ! level fully include in the ice shelf boundary layer 723 745 DO jk = ikt, ikb - 1 724 746 ze3 = fse3v_n(ji,jj,jk) 725 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3726 END DO 727 728 ! level partially include in ice shelf boundary layer747 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 748 END DO 749 750 ! level partially include in ice shelf boundary layer 729 751 zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 730 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,ikb) * (1._wp - zhk)752 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 731 753 END DO 732 754 END DO 733 755 DO jj = 2,jpj 734 756 DO ji = 2,jpi 735 varout(ji,jj) = 0.5_wp * (varout(ji,jj) + varout(ji,jj-1)) 736 END DO 737 END DO 738 CALL lbc_lnk(varout,'T',-1.) 739 END IF 740 741 ! compute T in the top boundary layer at T- point 742 IF (cptin == 'T') THEN 757 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 758 END DO 759 END DO 760 CALL lbc_lnk(pvarout,'T',-1.) 761 762 CASE ( 'T' ) ! compute T in the top boundary layer at T- point 743 763 DO jj = 1,jpj 744 764 DO ji = 1,jpi … … 746 766 ikb = misfkb(ji,jj) 747 767 748 ! level fully include in the ice shelf boundary layer768 ! level fully include in the ice shelf boundary layer 749 769 DO jk = ikt, ikb - 1 750 770 ze3 = fse3t_n(ji,jj,jk) 751 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3752 END DO 753 754 ! level partially include in ice shelf boundary layer771 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 772 END DO 773 774 ! level partially include in ice shelf boundary layer 755 775 zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 756 varout(ji,jj) = varout(ji,jj) +varin(ji,jj,ikb) * (1._wp - zhk)757 END DO 758 END DO 759 END IF776 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 777 END DO 778 END DO 779 END SELECT 760 780 761 781 ! mask mean tbl value 762 varout(:,:) =varout(:,:) * ssmask(:,:)782 pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 763 783 764 784 ! deallocation … … 780 800 !! ** Action : phdivn decreased by the runoff inflow 781 801 !!---------------------------------------------------------------------- 782 REAL(wp), DIMENSION(:,:,:), INTENT( inout) :: phdivn ! horizontal divergence783 ! !802 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdivn ! horizontal divergence 803 ! 784 804 INTEGER :: ji, jj, jk ! dummy loop indices 785 805 INTEGER :: ikt, ikb … … 790 810 zfact = 0.5_wp 791 811 ! 792 IF ( lk_vvl) THEN ! need to re compute level distribution of isf fresh water812 IF ( lk_vvl ) THEN ! need to re compute level distribution of isf fresh water 793 813 DO jj = 1,jpj 794 814 DO ji = 1,jpi … … 800 820 ! determine the deepest level influenced by the boundary layer 801 821 DO jk = ikt+1, mbkt(ji,jj) 802 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT.rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk822 IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 803 823 END DO 804 824 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 820 840 DO jk = ikt, ikb - 1 821 841 phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 822 & 842 & * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 823 843 END DO 824 844 ! level partially include in ice shelf boundary layer 825 845 phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 826 &+ fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)846 & + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 827 847 END DO 828 848 END DO
Note: See TracChangeset
for help on using the changeset viewer.
