Changeset 11217
- Timestamp:
- 2019-07-04T17:31:44+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
r11215 r11217 17 17 !! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/) 18 18 !! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 20 20 !!---------------------------------------------------------------------- 21 21 … … 24 24 !! sbc_blk : bulk formulation as ocean surface boundary condition 25 25 !! blk_oce : computes momentum, heat and freshwater fluxes over ocean 26 !! sea-ice case only : 26 !! sea-ice case only : 27 27 !! blk_ice_tau : provide the air-ice stress 28 28 !! blk_ice_flx : provide the heat and mass fluxes at air-ice interface 29 29 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 30 30 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 31 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 31 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 32 32 !!---------------------------------------------------------------------- 33 33 USE oce ! ocean dynamics and tracers … … 45 45 USE icethd_dh ! for CALL ice_thd_snwblow 46 46 #endif 47 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 48 USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003) 47 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 48 USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003) 49 49 USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013) 50 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31) 50 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31) 51 51 ! 52 52 USE iom ! I/O manager library … … 58 58 USE sbcblk_phy !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc... 59 59 60 60 61 61 IMPLICIT NONE 62 62 PRIVATE … … 68 68 PUBLIC blk_ice_flx ! routine called in icesbc 69 69 PUBLIC blk_ice_qcn ! routine called in icesbc 70 #endif 70 #endif 71 71 72 72 INTEGER , PARAMETER :: jpfld =10 ! maximum number of files to read … … 96 96 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 97 97 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 98 !!gm ref namelist initialize it so remove the setting to false below98 !!gm ref namelist initialize it so remove the setting to false below 99 99 LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 100 100 LOGICAL :: ln_Cd_L15 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) … … 111 111 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature if .true. (specific humidity espected if .false.) !LB 112 112 !LB. 113 113 114 114 INTEGER :: nblk ! choice of the bulk algorithm 115 115 ! ! associated indices: … … 159 159 !! ** Purpose : choose and initialize a bulk formulae formulation 160 160 !! 161 !! ** Method : 161 !! ** Method : 162 162 !! 163 163 !!---------------------------------------------------------------------- … … 173 173 & sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, & 174 174 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 175 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 175 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 176 176 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, & 177 177 & ln_skin, ln_humi_dpt ! cool-skin / warm-layer param. !LB … … 181 181 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 182 182 ! 183 ! !** read bulk namelist 183 ! !** read bulk namelist 184 184 REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters 185 185 READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) … … 194 194 ! !** initialization of the chosen bulk formulae (+ check) 195 195 ! !* select the bulk chosen in the namelist and check the choice 196 ioptio = 0 197 IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF 198 IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF 199 IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF 200 IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF 196 ioptio = 0 197 IF( ln_NCAR ) THEN 198 nblk = np_NCAR ; ioptio = ioptio + 1 199 ENDIF 200 IF( ln_COARE_3p0 ) THEN 201 nblk = np_COARE_3p0 ; ioptio = ioptio + 1 202 ENDIF 203 IF( ln_COARE_3p5 ) THEN 204 nblk = np_COARE_3p5 ; ioptio = ioptio + 1 205 ENDIF 206 IF( ln_ECMWF ) THEN 207 nblk = np_ECMWF ; ioptio = ioptio + 1 208 ENDIF 201 209 ! 202 210 !LB: … … 208 216 END IF 209 217 !LB. 210 218 211 219 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 212 220 ! 213 221 IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr 214 222 IF( sn_qsr%nfreqh /= 24 ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 215 IF( sn_qsr%ln_tint ) THEN 223 IF( sn_qsr%ln_tint ) THEN 216 224 CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', & 217 225 & ' ==> We force time interpolation = .false. for qsr' ) 218 sn_qsr%ln_tint = . false.226 sn_qsr%ln_tint = .FALSE. 219 227 ENDIF 220 228 ENDIF … … 245 253 ! 246 254 IF ( ln_wave ) THEN 247 !Activated wave module but neither drag nor stokes drift activated255 !Activated wave module but neither drag nor stokes drift activated 248 256 IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN 249 257 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 250 !drag coefficient read from wave model definable only with mfs bulk formulae and core251 ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR ) THEN 252 258 !drag coefficient read from wave model definable only with mfs bulk formulae and core 259 ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR ) THEN 260 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 253 261 ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN 254 262 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 255 263 ENDIF 256 264 ELSE 257 IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) &258 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', &259 & 'with drag coefficient (ln_cdgw =T) ' , &260 & 'or Stokes Drift (ln_sdw=T) ' , &261 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', &262 & 'or Stokes-Coriolis term (ln_stcori=T)' )263 ENDIF 264 ! 265 ! 265 IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 266 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 267 & 'with drag coefficient (ln_cdgw =T) ' , & 268 & 'or Stokes Drift (ln_sdw=T) ' , & 269 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', & 270 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 271 ENDIF 272 ! 273 ! 266 274 IF(lwp) THEN !** Control print 267 275 ! 268 WRITE(numout,*) !* namelist 276 WRITE(numout,*) !* namelist 269 277 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 270 278 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR … … 344 352 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 345 353 qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) 346 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 347 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 348 ENDIF 354 IF( ln_dm2dc ) THEN 355 qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 356 ELSE 357 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 358 ENDIF 349 359 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 350 360 !LB: … … 411 421 412 422 ! ... components ( U10m - U_oce ) at T-point (unmasked) 413 !!gm move zwnd_i (_j) set to zero inside the key_cyclone ???423 !!gm move zwnd_i (_j) set to zero inside the key_cyclone ??? 414 424 zwnd_i(:,:) = 0._wp 415 425 zwnd_j(:,:) = 0._wp … … 440 450 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 441 451 zztmp = 1. - albo 442 IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 443 ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 452 IF( ln_dm2dc ) THEN 453 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 454 ELSE 455 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 444 456 ENDIF 445 457 … … 451 463 ! ... specific humidity at SST and IST tmask( 452 464 zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 465 466 IF ( ln_skin ) THEN 467 zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action 468 zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action 469 END IF 453 470 454 471 !LB: … … 458 475 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !' 459 476 zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 460 ELSE 477 ELSE 461 478 zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 462 479 END IF … … 467 484 !! (since reanalysis products provide T at z, not theta !) 468 485 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt 469 470 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 471 ! 472 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! NCAR-COREv2 473 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 474 475 !LB: Skin!!! 476 CASE( np_COARE_3p0 ) 477 IF ( ln_skin ) THEN 486 487 488 IF ( ln_skin ) THEN 489 490 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 491 492 CASE( np_COARE_3p0 ) 478 493 CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.0 479 494 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 480 495 & Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 496 497 CASE( np_ECMWF ) 498 CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! ECMWF 499 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 500 & Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 501 502 CASE DEFAULT 503 CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' ) 504 END SELECT 505 506 !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and zsq: 507 WHERE ( fr_i < 0.001_wp ) 508 ! zst and zsq have been updated by cool-skin/warm-layer scheme and we keep it!!! 481 509 zst(:,:) = zst(:,:)*tmask(:,:,1) 482 510 zsq(:,:) = zsq(:,:)*tmask(:,:,1) 483 ELSE 484 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!' 511 ELSEWHERE 512 ! we forget about the update... 513 zst(:,:) = zqsb(:,:) !LB: using what we backed up before skin-algo 514 zsq(:,:) = zqla(:,:) !LB: " " " 515 END WHERE 516 517 !LB: Update of tsk, the official array for skin temperature 518 tsk(:,:) = zst(:,:) 519 520 ELSE !IF ( ln_skin ) 521 522 523 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 524 ! 525 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! NCAR-COREv2 526 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 527 528 CASE( np_COARE_3p0 ) 529 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare" WITHOUT CSWL optional arrays!!!' 485 530 CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.0 486 531 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 487 END IF 488 489 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.5 490 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 491 492 !LB: Skin!!! 493 CASE( np_ECMWF ) 494 IF ( ln_skin ) THEN 495 CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! ECMWF 496 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 497 & Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 498 zst(:,:) = zst(:,:)*tmask(:,:,1) 499 zsq(:,:) = zsq(:,:)*tmask(:,:,1) 500 ELSE 532 533 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.5 534 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 535 536 CASE( np_ECMWF ) 501 537 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' 502 CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! ECMWF 538 CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! ECMWF 503 539 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 504 END IF 505 !LB. 506 507 CASE DEFAULT 508 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 509 END SELECT 510 511 512 !LB: cool-skin/warm-layer: 513 IF ( ln_skin ) tsk(:,:) = zst(:,:) 540 541 CASE DEFAULT 542 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 543 END SELECT 544 545 END IF ! IF ( ln_skin ) 546 547 514 548 515 549 … … 521 555 END IF 522 556 523 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef.524 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef.557 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. 558 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 525 559 526 560 DO jj = 1, jpj ! tau module, i and j component … … 612 646 ! 613 647 !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 648 CALL iom_put( "rho_air" , rhoa ) ! output air density (kg/m^3) !#LB 649 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 650 CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean 651 CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean 652 CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean 653 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 654 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 655 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 656 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 657 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 658 CALL iom_put( 'snowpre', sprecip ) ! Snow 659 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 660 IF ( ln_skin ) THEN 661 CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius 662 CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference... 663 END IF 630 664 !!#LB. ENDIF 631 665 ! … … 650 684 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 651 685 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 652 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 686 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 653 687 !!---------------------------------------------------------------------- 654 688 … … 693 727 Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical 694 728 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 695 CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 696 ENDIF 697 698 !! CALL iom_put( "rCd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef.699 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef.729 CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 730 ENDIF 731 732 !! CALL iom_put( "rCd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef. 733 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef. 700 734 701 735 ! local scalars ( place there for vector optimisation purposes) … … 776 810 777 811 zztmp = 1. / ( 1. - albo ) 778 WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 779 ELSEWHERE ; z1_st(:,:,:) = 0._wp 812 WHERE( ptsu(:,:,:) /= 0._wp ) 813 z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 814 ELSEWHERE 815 z1_st(:,:,:) = 0._wp 780 816 END WHERE 781 817 ! ! ========================== ! … … 866 902 DO jl = 1, jpl 867 903 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 868 ! ! But we do not have Tice => consider it at 0degC => evap=0 904 ! ! But we do not have Tice => consider it at 0degC => evap=0 869 905 END DO 870 906 … … 873 909 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 874 910 ! 875 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 911 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 876 912 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 877 913 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 878 914 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 879 915 ELSEWHERE ! zero when hs>0 880 qtr_ice_top(:,:,:) = 0._wp 916 qtr_ice_top(:,:,:) = 0._wp 881 917 END WHERE 882 918 ! … … 891 927 ! 892 928 END SUBROUTINE blk_ice_flx 893 929 894 930 895 931 SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) … … 900 936 !! to force sea ice / snow thermodynamics 901 937 !! in the case conduction flux is emulated 902 !! 938 !! 903 939 !! ** Method : compute surface energy balance assuming neglecting heat storage 904 940 !! following the 0-layer Semtner (1976) approach … … 925 961 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor 926 962 !!--------------------------------------------------------------------- 927 963 928 964 ! -------------------------------------! 929 965 ! I Enhanced conduction factor ! … … 933 969 ! 934 970 zgfac(:,:,:) = 1._wp 935 971 936 972 IF( ld_virtual_itd ) THEN 937 973 ! … … 939 975 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 940 976 zfac3 = 2._wp / zepsilon 941 ! 942 DO jl = 1, jpl 977 ! 978 DO jl = 1, jpl 943 979 DO jj = 1 , jpj 944 980 DO ji = 1, jpi … … 948 984 END DO 949 985 END DO 950 ! 951 ENDIF 952 986 ! 987 ENDIF 988 953 989 ! -------------------------------------------------------------! 954 990 ! II Surface temperature and conduction flux ! … … 960 996 DO jj = 1 , jpj 961 997 DO ji = 1, jpi 962 ! 998 ! 963 999 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 964 1000 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) … … 977 1013 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 978 1014 qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) ) & 979 1015 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 980 1016 981 1017 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 982 hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 1018 hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 983 1019 984 1020 END DO 985 1021 END DO 986 1022 ! 987 END DO 988 ! 1023 END DO 1024 ! 989 1025 END SUBROUTINE blk_ice_qcn 990 1026 991 1027 992 1028 SUBROUTINE Cdn10_Lupkes2012( Cd ) … … 994 1030 !! *** ROUTINE Cdn10_Lupkes2012 *** 995 1031 !! 996 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1032 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 997 1033 !! to make it dependent on edges at leads, melt ponds and flows. 998 1034 !! After some approximations, this can be resumed to a dependency 999 1035 !! on ice concentration. 1000 !! 1036 !! 1001 1037 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50) 1002 1038 !! with the highest level of approximation: level4, eq.(59) … … 1010 1046 !! 1011 1047 !! This new drag has a parabolic shape (as a function of A) starting at 1012 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1048 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1013 1049 !! and going down to Cdi(say 1.4e-3) for A=1 1014 1050 !! … … 1037 1073 Cd(:,:) = rCd_ice + & ! pure ice drag 1038 1074 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1039 1075 1040 1076 END SUBROUTINE Cdn10_Lupkes2012 1041 1077 … … 1046 1082 !! 1047 1083 !! ** pUrpose : Alternative turbulent transfert coefficients formulation 1048 !! between sea-ice and atmosphere with distinct momentum 1049 !! and heat coefficients depending on sea-ice concentration 1084 !! between sea-ice and atmosphere with distinct momentum 1085 !! and heat coefficients depending on sea-ice concentration 1050 1086 !! and atmospheric stability (no meltponds effect for now). 1051 !! 1087 !! 1052 1088 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1053 1089 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1054 1090 !! it considers specific skin and form drags (Andreas et al. 2010) 1055 !! to compute neutral transfert coefficients for both heat and 1091 !! to compute neutral transfert coefficients for both heat and 1056 1092 !! momemtum fluxes. Atmospheric stability effect on transfert 1057 1093 !! coefficient is also taken into account following Louis (1979). … … 1093 1129 1094 1130 ! mean temperature 1095 WHERE( at_i_b(:,:) > 1.e-20 ) ; ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:) 1096 ELSEWHERE ; ztm_su(:,:) = rt0 1131 WHERE( at_i_b(:,:) > 1.e-20 ) 1132 ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:) 1133 ELSEWHERE 1134 ztm_su(:,:) = rt0 1097 1135 ENDWHERE 1098 1136 1099 1137 ! Momentum Neutral Transfert Coefficients (should be a constant) 1100 1138 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 … … 1105 1143 ! Heat Neutral Transfert Coefficients 1106 1144 zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) ) ! Eq. 50 + Eq. 52 (cf Lupkes email for details) 1107 1145 1108 1146 ! Atmospheric and Surface Variables 1109 1147 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin … … 1117 1155 zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1118 1156 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1119 1157 1120 1158 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1121 1159 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1122 1160 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1123 1161 1124 1162 ! Momentum and Heat Neutral Transfert Coefficients 1125 1163 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1126 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1127 1164 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1165 1128 1166 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 1129 1167 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water … … 1137 1175 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1138 1176 ENDIF 1139 1177 1140 1178 IF( zrib_i <= 0._wp ) THEN 1141 1179 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 … … 1145 1183 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1146 1184 ENDIF 1147 1185 1148 1186 ! Momentum Transfert Coefficients (Eq. 38) 1149 1187 Cd(ji,jj) = zCdn_skin_ice * zfmi + & 1150 1188 & zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 1151 1189 1152 1190 ! Heat Transfert Coefficients (Eq. 49) 1153 1191 Ch(ji,jj) = zChn_skin_ice * zfhi + &
Note: See TracChangeset
for help on using the changeset viewer.
