- Timestamp:
- 2017-12-13T15:58:53+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
r8329 r9019 5 5 !! shelf 6 6 !!====================================================================== 7 !! History : 3.2 8 !! X.X 9 !! 3.4 7 !! History : 3.2 ! 2011-02 (C.Harris ) Original code isf cav 8 !! X.X ! 2006-02 (C. Wang ) Original code bg03 9 !! 3.4 ! 2013-03 (P. Mathiot) Merging + parametrization 10 10 !!---------------------------------------------------------------------- 11 11 12 12 !!---------------------------------------------------------------------- 13 !! sbc_isf 13 !! sbc_isf : update sbc under ice shelf 14 14 !!---------------------------------------------------------------------- 15 USE oce 16 USE dom_oce 17 USE phycst 18 USE eosbn2 19 USE sbc_oce 20 USE zdf bfr !15 USE oce ! ocean dynamics and tracers 16 USE dom_oce ! ocean space and time domain 17 USE phycst ! physical constants 18 USE eosbn2 ! equation of state 19 USE sbc_oce ! surface boundary condition: ocean fields 20 USE zdfdrg ! vertical physics: top/bottom drag coef. 21 21 ! 22 USE in_out_manager ! I/O manager 23 USE iom ! I/O manager library 24 USE fldread ! read input field at current time step 25 USE lbclnk ! 26 USE wrk_nemo ! Memory allocation 27 USE timing ! Timing 28 USE lib_fortran ! glob_sum 22 USE in_out_manager ! I/O manager 23 USE iom ! I/O library 24 USE fldread ! read input field at current time step 25 USE lbclnk ! 26 USE timing ! Timing 27 USE lib_fortran ! glob_sum 29 28 30 29 IMPLICIT NONE … … 35 34 ! public in order to be able to output then 36 35 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc !: before and now T & S isf contents [K.m/s & PSU.m/s]38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf [W/m2]39 36 REAL(wp), PUBLIC :: rn_hisf_tbl !: thickness of top boundary layer [m] 40 37 INTEGER , PUBLIC :: nn_isf !: flag to choose between explicit/param/specified … … 44 41 REAL(wp), PUBLIC :: rn_gammas0 !: salinity exchange coeficient [] 45 42 46 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rzisf_tbl !:depth of calving front (shallowest point) nn_isf ==2/3 47 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: rhisf_tbl, rhisf_tbl_0 !:thickness of tbl [m] 48 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: r1_hisf_tbl !:1/thickness of tbl 49 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ralpha !:proportion of bottom cell influenced by tbl 50 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: risfLeff !:effective length (Leff) BG03 nn_isf==2 51 REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 52 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:) :: misfkt, misfkb !:Level of ice shelf base 53 54 LOGICAL, PUBLIC :: l_isfcpl = .false. ! isf recieved from oasis 55 56 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp ! specific heat of ice shelf [J/kg/K] 57 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp ! heat diffusivity through the ice-shelf [m2/s] 58 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp ! volumic mass of ice shelf [kg/m3] 59 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp ! air temperature on top of ice shelf [C] 60 REAL(wp), PUBLIC, SAVE :: rlfusisf = 0.334e6_wp ! latent heat of fusion of ice shelf [J/kg] 43 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt , misfkb !: Level of ice shelf base 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rzisf_tbl !: depth of calving front (shallowest point) nn_isf ==2/3 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl, rhisf_tbl_0 !: thickness of tbl [m] 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hisf_tbl !: 1/thickness of tbl 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ralpha !: proportion of bottom cell influenced by tbl 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfLeff !: effective length (Leff) BG03 nn_isf==2 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ttbl, stbl, utbl, vtbl !: top boundary layer variable at T point 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: qisf !: net heat flux from ice shelf [W/m2] 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_tsc_b, risf_tsc !: before and now T & S isf contents [K.m/s & PSU.m/s] 52 53 LOGICAL, PUBLIC :: l_isfcpl = .false. !: isf recieved from oasis 54 55 REAL(wp), PUBLIC, SAVE :: rcpi = 2000.0_wp !: specific heat of ice shelf [J/kg/K] 56 REAL(wp), PUBLIC, SAVE :: rkappa = 1.54e-6_wp !: heat diffusivity through the ice-shelf [m2/s] 57 REAL(wp), PUBLIC, SAVE :: rhoisf = 920.0_wp !: volumic mass of ice shelf [kg/m3] 58 REAL(wp), PUBLIC, SAVE :: tsurf = -20.0_wp !: air temperature on top of ice shelf [C] 59 REAL(wp), PUBLIC, SAVE :: rlfusisf = 0.334e6_wp !: latent heat of fusion of ice shelf [J/kg] 61 60 62 61 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) … … 71 70 72 71 !!---------------------------------------------------------------------- 73 !! NEMO/OPA 3.7 , LOCEAN-IPSL (2015)72 !! NEMO/OPA 4.0 , LOCEAN-IPSL (2017) 74 73 !! $Id$ 75 74 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 77 76 CONTAINS 78 77 79 SUBROUTINE sbc_isf( kt)78 SUBROUTINE sbc_isf( kt ) 80 79 !!--------------------------------------------------------------------- 81 80 !! *** ROUTINE sbc_isf *** … … 90 89 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 91 90 !!---------------------------------------------------------------------- 92 INTEGER, INTENT( in ) :: kt! ocean time step93 ! 94 INTEGER :: ji, jj, jk! loop index95 INTEGER :: ikt, ikb ! loop index96 REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep! freezing temperature (zt_frz) at depth (zdep)97 REAL(wp), DIMENSION(:,: ,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d98 REAL(wp), DIMENSION(:,: ), POINTER :: zqhcisf2d91 INTEGER, INTENT(in) :: kt ! ocean time step 92 ! 93 INTEGER :: ji, jj, jk ! loop index 94 INTEGER :: ikt, ikb ! local integers 95 REAL(wp), DIMENSION(jpi,jpj) :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 96 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zqhcisf2d 97 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zfwfisf3d, zqhcisf3d, zqlatisf3d 99 98 !!--------------------------------------------------------------------- 100 99 ! 101 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 102 ! allocation 103 CALL wrk_alloc( jpi,jpj, zt_frz, zdep ) 104 105 ! compute salt and heat flux 100 IF( MOD( kt-1, nn_fsbc) == 0 ) THEN ! compute salt and heat flux 101 ! 106 102 SELECT CASE ( nn_isf ) 107 103 CASE ( 1 ) ! realistic ice shelf formulation … … 121 117 ELSE ; qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 122 118 ENDIF 123 119 ! 124 120 CASE ( 2 ) ! Beckmann and Goosse parametrisation 125 121 stbl(:,:) = soce 126 122 CALL sbc_isf_bg03(kt) 127 123 ! 128 124 CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation) 129 125 ! specified runoff in depth (Mathiot et al., XXXX in preparation) … … 134 130 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 135 131 stbl(:,:) = soce 136 132 ! 137 133 CASE ( 4 ) ! specified fwf and heat flux forcing beneath the ice shelf 138 ! specified fwf and heat flux forcing beneath the ice shelf134 ! ! specified fwf and heat flux forcing beneath the ice shelf 139 135 IF( .NOT.l_isfcpl ) THEN 140 136 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) … … 144 140 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 145 141 stbl(:,:) = soce 146 142 ! 147 143 END SELECT 148 144 … … 162 158 163 159 ! lbclnk 164 CALL lbc_lnk( risf_tsc(:,:,jp_tem),'T',1.)165 CALL lbc_lnk( risf_tsc(:,:,jp_sal),'T',1.)166 CALL lbc_lnk( fwfisf(:,:),'T',1.)167 CALL lbc_lnk( qisf(:,:),'T',1.)160 CALL lbc_lnk( risf_tsc(:,:,jp_tem),'T',1.) 161 CALL lbc_lnk( risf_tsc(:,:,jp_sal),'T',1.) 162 CALL lbc_lnk( fwfisf (:,:) ,'T',1.) 163 CALL lbc_lnk( qisf (:,:) ,'T',1.) 168 164 169 165 ! output … … 173 169 IF( iom_use('fwfisf' ) ) CALL iom_put( 'fwfisf' , fwfisf(:,:) ) ! isf mass flux (opposite sign) 174 170 175 ! Diagnostics176 IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN177 CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d)178 CALL wrk_alloc( jpi,jpj, zqhcisf2d)179 180 zfwfisf3d (:,:,:) = 0.0_wp ! 3d ice shelf melting (kg/m2/s)181 zqhcisf3d (:,:,:) = 0.0_wp ! 3d heat content flux (W/m2)182 zqlatisf3d(:,:,:) = 0.0_wp ! 3d ice shelf melting latent heat flux (W/m2)183 zqhcisf2d (:,:) = fwfisf(:,:) * zt_frz * rcp! 2d heat content flux (W/m2)184 171 ! Diagnostics 172 IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN 173 ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) ) 174 ALLOCATE( zqhcisf2d(jpi,jpj) ) 175 ! 176 zfwfisf3d (:,:,:) = 0._wp ! 3d ice shelf melting (kg/m2/s) 177 zqhcisf3d (:,:,:) = 0._wp ! 3d heat content flux (W/m2) 178 zqlatisf3d(:,:,:) = 0._wp ! 3d ice shelf melting latent heat flux (W/m2) 179 zqhcisf2d (:,:) = fwfisf(:,:) * zt_frz * rcp ! 2d heat content flux (W/m2) 180 ! 185 181 DO jj = 1,jpj 186 182 DO ji = 1,jpi … … 200 196 END DO 201 197 END DO 202 198 ! 203 199 CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:)) 204 200 CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:)) 205 201 CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:)) 206 202 CALL iom_put('qhcisf' , zqhcisf2d (:,: )) 207 208 CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d ) 209 CALL wrk_dealloc( jpi,jpj, zqhcisf2d ) 210 END IF 211 ! deallocation 212 CALL wrk_dealloc( jpi,jpj, zt_frz, zdep ) 213 ! 214 END IF 215 216 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 217 IF( ln_rstart .AND. & ! Restart: read in restart file 218 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 219 IF(lwp) WRITE(numout,*) ' nit000-1 isf tracer content forcing fields read in the restart file' 220 CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) ) ! before salt content isf_tsc trend 221 CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) ) ! before salt content isf_tsc trend 222 CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) ) ! before salt content isf_tsc trend 223 ELSE 224 fwfisf_b(:,:) = fwfisf(:,:) 225 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 226 END IF 227 END IF 228 ! 229 IF( lrst_oce ) THEN 230 IF(lwp) WRITE(numout,*) 231 IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ', & 232 & 'at it= ', kt,' date= ', ndastp 233 IF(lwp) WRITE(numout,*) '~~~~' 234 CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) ) 235 CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) ) 236 CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) ) 203 ! 204 DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d ) 205 DEALLOCATE( zqhcisf2d ) 237 206 ENDIF 238 207 ! 239 END SUBROUTINE sbc_isf 240 241 242 INTEGER FUNCTION sbc_isf_alloc() 208 ENDIF 209 210 IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! 211 IF( ln_rstart .AND. & ! Restart: read in restart file 212 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 213 IF(lwp) WRITE(numout,*) ' nit000-1 isf tracer content forcing fields read in the restart file' 214 CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) ) ! before salt content isf_tsc trend 215 CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) ) ! before salt content isf_tsc trend 216 CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) ) ! before salt content isf_tsc trend 217 ELSE 218 fwfisf_b(:,:) = fwfisf(:,:) 219 risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 220 ENDIF 221 ENDIF 222 ! 223 IF( lrst_oce ) THEN 224 IF(lwp) WRITE(numout,*) 225 IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ', & 226 & 'at it= ', kt,' date= ', ndastp 227 IF(lwp) WRITE(numout,*) '~~~~' 228 CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) ) 229 CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) ) 230 CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) ) 231 ENDIF 232 ! 233 END SUBROUTINE sbc_isf 234 235 236 INTEGER FUNCTION sbc_isf_alloc() 243 237 !!---------------------------------------------------------------------- 244 238 !! *** FUNCTION sbc_isf_rnf_alloc *** … … 256 250 IF( sbc_isf_alloc /= 0 ) CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 257 251 ! 258 END IF 259 END FUNCTION 252 ENDIF 253 END FUNCTION 254 260 255 261 256 SUBROUTINE sbc_isf_init … … 294 289 295 290 IF ( lwp ) WRITE(numout,*) 296 IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 297 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 298 IF ( lwp ) WRITE(numout,*) 'sbcisf :' 299 IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 291 IF ( lwp ) WRITE(numout,*) 'sbc_isf_init : heat flux of the ice shelf' 292 IF ( lwp ) WRITE(numout,*) '~~~~~~~~~~~' 300 293 IF ( lwp ) WRITE(numout,*) ' nn_isf = ', nn_isf 301 294 IF ( lwp ) WRITE(numout,*) ' nn_isfblk = ', nn_isfblk … … 304 297 IF ( lwp ) WRITE(numout,*) ' rn_gammat0 = ', rn_gammat0 305 298 IF ( lwp ) WRITE(numout,*) ' rn_gammas0 = ', rn_gammas0 306 IF ( lwp ) WRITE(numout,*) ' rn_ tfri2 = ', rn_tfri2299 IF ( lwp ) WRITE(numout,*) ' rn_Cd0 = ', r_Cdmin_top 307 300 ! 308 301 ! Allocate public variable … … 310 303 ! 311 304 ! initialisation 312 qisf (:,:) = 0._wp ;fwfisf (:,:) = 0._wp313 risf_tsc(:,:,:) = 0._wp ;fwfisf_b(:,:) = 0._wp305 qisf (:,:) = 0._wp ; fwfisf (:,:) = 0._wp 306 risf_tsc(:,:,:) = 0._wp ; fwfisf_b(:,:) = 0._wp 314 307 ! 315 308 ! define isf tbl tickness, top and bottom indice … … 317 310 CASE ( 1 ) 318 311 rhisf_tbl(:,:) = rn_hisf_tbl 319 misfkt (:,:)= mikt(:,:) ! same indice for bg03 et cav => used in isfdiv312 misfkt (:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 320 313 321 314 CASE ( 2 , 3 ) … … 351 344 DO jj = 1, jpj 352 345 ik = 2 346 !!gm potential bug: use gdepw_0 not _n 353 347 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO 354 348 misfkt(ji,jj) = ik-1 … … 359 353 ! as in nn_isf == 1 360 354 rhisf_tbl(:,:) = rn_hisf_tbl 361 misfkt (:,:)= mikt(:,:) ! same indice for bg03 et cav => used in isfdiv355 misfkt (:,:) = mikt(:,:) ! same indice for bg03 et cav => used in isfdiv 362 356 363 357 ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) … … 382 376 ! determine the deepest level influenced by the boundary layer 383 377 DO jk = ikt+1, mbkt(ji,jj) 384 IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )ikb = jk378 IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 385 379 END DO 386 380 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 395 389 END SUBROUTINE sbc_isf_init 396 390 391 397 392 SUBROUTINE sbc_isf_bg03(kt) 398 393 !!--------------------------------------------------------------------- … … 407 402 !! interaction for climate models", Ocean Modelling 5(2003) 157-170. 408 403 !! (hereafter BG) 409 !! History : 410 !! 06-02 (C. Wang) Original code 404 !! History : 06-02 (C. Wang) Original code 411 405 !!---------------------------------------------------------------------- 412 406 INTEGER, INTENT ( in ) :: kt … … 420 414 !!---------------------------------------------------------------------- 421 415 422 IF ( nn_timing == 1 )CALL timing_start('sbc_isf_bg03')416 IF( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 423 417 ! 424 418 DO ji = 1, jpi … … 446 440 !add to salinity trend 447 441 ELSE 448 qisf(ji,jj) = 0._wp ;fwfisf(ji,jj) = 0._wp442 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 449 443 END IF 450 444 END DO 451 445 END DO 452 446 ! 453 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03')447 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_bg03') 454 448 ! 455 449 END SUBROUTINE sbc_isf_bg03 450 456 451 457 452 SUBROUTINE sbc_isf_cav( kt ) … … 468 463 !! emp, emps : update freshwater flux below ice shelf 469 464 !!--------------------------------------------------------------------- 470 INTEGER, INTENT(in) :: kt! ocean time step465 INTEGER, INTENT(in) :: kt ! ocean time step 471 466 ! 472 467 INTEGER :: ji, jj ! dummy loop indices 473 468 INTEGER :: nit 469 LOGICAL :: lit 474 470 REAL(wp) :: zlamb1, zlamb2, zlamb3 475 471 REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 … … 477 473 REAL(wp) :: zeps = 1.e-20_wp 478 474 REAL(wp) :: zerr 479 REAL(wp), DIMENSION(:,:), POINTER :: zfrz 480 REAL(wp), DIMENSION(:,:), POINTER :: zgammat, zgammas 481 REAL(wp), DIMENSION(:,:), POINTER :: zfwflx, zhtflx, zhtflx_b 482 LOGICAL :: lit 475 REAL(wp), DIMENSION(jpi,jpj) :: zfrz 476 REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas 477 REAL(wp), DIMENSION(jpi,jpj) :: zfwflx, zhtflx, zhtflx_b 483 478 !!--------------------------------------------------------------------- 484 479 ! coeficient for linearisation of potential tfreez … … 489 484 IF( nn_timing == 1 ) CALL timing_start('sbc_isf_cav') 490 485 ! 491 CALL wrk_alloc( jpi,jpj, zfrz , zgammat, zgammas )492 CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )493 494 486 ! initialisation 495 487 zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 … … 583 575 CALL iom_put('isfgammas', zgammas) 584 576 ! 585 CALL wrk_dealloc( jpi,jpj, zfrz , zgammat, zgammas )586 CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )587 !588 577 IF( nn_timing == 1 ) CALL timing_stop('sbc_isf_cav') 589 578 ! … … 605 594 INTEGER :: ikt 606 595 INTEGER :: ji, jj ! loop index 607 REAL(wp), DIMENSION(:,:), POINTER :: zustar ! U, V at T point and friction velocity608 596 REAL(wp) :: zdku, zdkv ! U, V shear 609 597 REAL(wp) :: zPr, zSc, zRc ! Prandtl, Scmidth and Richardson number … … 619 607 REAL(wp), PARAMETER :: znu = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 620 608 REAL(wp), DIMENSION(2) :: zts, zab 609 REAL(wp), DIMENSION(jpi,jpj) :: zustar ! U, V at T point and friction velocity 621 610 !!--------------------------------------------------------------------- 622 CALL wrk_alloc( jpi,jpj, zustar )623 611 ! 624 612 SELECT CASE ( nn_gammablk ) … … 631 619 !! Jenkins et al., 2010, JPO, p2298-2312 632 620 !! Adopted by Asay-Davis et al. (2015) 633 634 !! compute ustar (eq. 24) 635 zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 621 !!gm I don't understand the u* expression in those papers... (see for example zdfglf module) 622 !! for me ustar= Cd0 * |U| not (Cd0)^1/2 * |U| .... which is what you can find in Jenkins et al. 623 624 !! compute ustar (eq. 24) !! NB: here r_Cdmin_top = rn_Cd0 read in namdrg_top namelist) 625 zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 636 626 637 627 !! Compute gammats … … 643 633 !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 644 634 !! compute ustar 645 zustar(:,:) = SQRT( r n_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )635 zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) ) 646 636 647 637 !! compute Pr and Sc number (can be improved) … … 654 644 655 645 !! compute gamma 656 DO ji =2,jpi657 DO jj =2,jpj646 DO ji = 2, jpi 647 DO jj = 2, jpj 658 648 ikt = mikt(ji,jj) 659 649 660 IF (zustar(ji,jj) == 0._wp) THEN ! only for kt = 1 I think650 IF( zustar(ji,jj) == 0._wp ) THEN ! only for kt = 1 I think 661 651 pgt = rn_gammat0 662 652 pgs = rn_gammas0 663 653 ELSE 664 654 !! compute Rc number (as done in zdfric.F90) 655 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation 656 !!gm moreover, use Max(rn2,0) to take care of static instabilities.... 665 657 zcoef = 0.5_wp / e3w_n(ji,jj,ikt) 666 658 ! ! shear of horizontal velocity … … 708 700 CALL lbc_lnk(pgs(:,:),'T',1.) 709 701 END SELECT 710 CALL wrk_dealloc( jpi,jpj, zustar )711 702 ! 712 703 END SUBROUTINE sbc_isf_gammats 713 704 705 714 706 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 715 707 !!---------------------------------------------------------------------- … … 719 711 !! 720 712 !!---------------------------------------------------------------------- 721 REAL(wp), DIMENSION(:,:,:), INTENT( in ) :: pvarin 722 REAL(wp), DIMENSION(:,:) , INTENT( out ) :: pvarout 723 CHARACTER(len=1), INTENT( in ) :: cd_ptin ! point of variable in/out 724 ! 725 REAL(wp) :: ze3, zhk 726 REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 727 728 INTEGER :: ji, jj, jk ! loop index 729 INTEGER :: ikt, ikb ! top and bottom index of the tbl 730 !!---------------------------------------------------------------------- 731 ! allocation 732 CALL wrk_alloc( jpi,jpj, zhisf_tbl) 713 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pvarin 714 REAL(wp), DIMENSION(:,:) , INTENT( out) :: pvarout 715 CHARACTER(len=1), INTENT(in ) :: cd_ptin ! point of variable in/out 716 ! 717 INTEGER :: ji, jj, jk ! loop index 718 INTEGER :: ikt, ikb ! top and bottom index of the tbl 719 REAL(wp) :: ze3, zhk 720 REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl 721 !!---------------------------------------------------------------------- 733 722 734 723 ! initialisation … … 741 730 ikt = miku(ji,jj) ; ikb = miku(ji,jj) 742 731 ! thickness of boundary layer at least the top level thickness 743 zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj), e3u_n(ji,jj,ikt))732 zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) ) 744 733 745 734 ! determine the deepest level influenced by the boundary layer … … 760 749 END DO 761 750 END DO 762 DO jj = 2,jpj 763 DO ji = 2,jpi 751 DO jj = 2, jpj 752 DO ji = 2, jpi 753 !!gm a wet-point only average should be used here !!! 764 754 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 765 755 END DO … … 791 781 END DO 792 782 END DO 793 DO jj = 2,jpj 794 DO ji = 2,jpi 783 DO jj = 2, jpj 784 DO ji = 2, jpi 785 !!gm a wet-point only average should be used here !!! 795 786 pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 796 787 END DO … … 816 807 END DO 817 808 END SELECT 818 809 ! 819 810 ! mask mean tbl value 820 811 pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 821 822 ! deallocation823 CALL wrk_dealloc( jpi,jpj, zhisf_tbl )824 812 ! 825 813 END SUBROUTINE sbc_isf_tbl … … 887 875 ! 888 876 END SUBROUTINE sbc_isf_div 877 889 878 !!====================================================================== 890 879 END MODULE sbcisf
Note: See TracChangeset
for help on using the changeset viewer.
