Changeset 13655 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90
- Timestamp:
- 2020-10-21T16:15:13+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90
r13501 r13655 30 30 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 31 31 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 32 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag33 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag34 32 !!---------------------------------------------------------------------- 35 33 USE oce ! ocean dynamics and tracers … … 41 39 USE sbcdcy ! surface boundary condition: diurnal cycle 42 40 USE sbcwave , ONLY : cdn_wave ! wave module 43 USE sbc_ice ! Surface boundary condition: ice fields44 41 USE lib_fortran ! to use key_nosignedzero 42 ! 45 43 #if defined key_si3 44 USE sbc_ice ! Surface boundary condition: ice fields #LB? ok to be in 'key_si3' ??? 46 45 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif, nn_qtrice 47 46 USE icevar ! for CALL ice_var_snwblow 48 #endif 49 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 47 USE sbcblk_algo_ice_lu12 48 USE sbcblk_algo_ice_lg15 49 #endif 50 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - (formerly known as CORE, Large & Yeager, 2009) 50 51 USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 51 52 USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 52 53 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 45r1) 54 USE sbcblk_algo_andreas ! => turb_andreas : Andreas et al. 2015 55 ! 56 53 57 ! 54 58 USE iom ! I/O manager library … … 58 62 USE prtctl ! Print control 59 63 60 USE sbc blk_phy! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc...64 USE sbc_phy ! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc... 61 65 62 66 … … 100 104 LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013) 101 105 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1) 106 LOGICAL :: ln_ANDREAS ! "ANDREAS" algorithm (Andreas et al. 2015) 102 107 ! 103 LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 104 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 108 !#LB: 109 LOGICAL :: ln_Cx_ice_cst ! use constant ice-air bulk transfer coefficients (value given in namelist's rn_Cd_i, rn_Ce_i & rn_Ch_i) 110 REAL(wp) :: rn_Cd_i, rn_Ce_i, rn_Ch_i 111 LOGICAL :: ln_Cx_ice_LU12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 112 LOGICAL :: ln_Cx_ice_LG15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 113 !#LB. 105 114 ! 106 115 LOGICAL :: ln_crt_fbk ! Add surface current feedback to the wind stress computation (Renault et al. 2020) 107 116 REAL(wp) :: rn_stau_a ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta 108 REAL(wp) :: rn_stau_b ! 117 REAL(wp) :: rn_stau_b ! 109 118 ! 110 119 REAL(wp) :: rn_pfac ! multiplication factor for precipitation … … 113 122 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 114 123 ! 115 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme and ABL) 116 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice 117 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 124 INTEGER :: nn_iter_algo ! Number of iterations in bulk param. algo ("stable ABL + weak wind" requires more) 125 126 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: theta_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 127 128 #if defined key_si3 129 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice !#LB transfert coefficients over ice 130 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: theta_zu_i, q_zu_i !#LB fixme ! air temp. and spec. hum. over ice at wind speed height (L15 bulk scheme) 131 #endif 132 118 133 119 134 LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB … … 136 151 INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013) 137 152 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 45r1) 153 INTEGER, PARAMETER :: np_ANDREAS = 5 ! "ANDREAS" algorithm (Andreas et al. 2015) 154 155 !#LB: 156 #if defined key_si3 157 ! Same, over sea-ice: 158 INTEGER :: nblk_ice ! choice of the bulk algorithm 159 ! ! associated indices: 160 INTEGER, PARAMETER :: np_ice_cst = 1 ! constant transfer coefficients 161 INTEGER, PARAMETER :: np_ice_lu12 = 2 ! Lupkes, 2012 162 INTEGER, PARAMETER :: np_ice_lg15 = 3 ! Lupkes & Gryanik, 2015 163 #endif 164 !LB. 165 166 138 167 139 168 !! * Substitutions … … 150 179 !! *** ROUTINE sbc_blk_alloc *** 151 180 !!------------------------------------------------------------------- 152 ALLOCATE( t_zu(jpi,jpj) , q_zu(jpi,jpj) , & 153 & Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj), & 154 & Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc ) 155 ! 181 ALLOCATE( theta_zu(jpi,jpj), q_zu(jpi,jpj), STAT=sbc_blk_alloc ) 156 182 CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) 157 183 IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' ) 158 184 END FUNCTION sbc_blk_alloc 185 186 #if defined key_si3 187 INTEGER FUNCTION sbc_blk_ice_alloc() 188 !!------------------------------------------------------------------- 189 !! *** ROUTINE sbc_blk_ice_alloc *** 190 !!------------------------------------------------------------------- 191 ALLOCATE( Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), & 192 & theta_zu_i(jpi,jpj), q_zu_i(jpi,jpj), STAT=sbc_blk_ice_alloc ) 193 CALL mpp_sum ( 'sbcblk', sbc_blk_ice_alloc ) 194 IF( sbc_blk_ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_ice_alloc: failed to allocate arrays' ) 195 END FUNCTION sbc_blk_ice_alloc 196 #endif 159 197 160 198 … … 178 216 TYPE(FLD_N) :: sn_cc, sn_hpgi, sn_hpgj ! " " 179 217 INTEGER :: ipka ! number of levels in the atmospheric variable 180 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 181 & sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm, & 182 & sn_cc, sn_hpgi, sn_hpgj, & 183 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 184 & cn_dir , rn_zqt, rn_zu, & 185 & rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15, ln_tpot, & 186 & ln_crt_fbk, rn_stau_a, rn_stau_b, & ! current feedback 187 & ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB 218 NAMELIST/namsbc_blk/ ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, ln_ANDREAS, & ! bulk algorithm 219 & rn_zqt, rn_zu, nn_iter_algo, ln_skin_cs, ln_skin_wl, & 220 & rn_pfac, rn_efac, & 221 & ln_crt_fbk, rn_stau_a, rn_stau_b, & ! current feedback 222 & ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tpot, & 223 & ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i, & 224 & ln_Cx_ice_LU12, ln_Cx_ice_LG15, & 225 & cn_dir, & 226 & sn_wndi, sn_wndj, sn_qsr, sn_qlw , & ! input fields 227 & sn_tair, sn_humi, sn_prec, sn_snow, sn_slp, & 228 & sn_uoatm, sn_voatm, sn_cc, sn_hpgi, sn_hpgj 229 230 ! cool-skin / warm-layer !LB 188 231 !!--------------------------------------------------------------------- 189 232 ! 190 233 ! ! allocate sbc_blk_core array 191 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 234 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 235 ! 236 #if defined key_si3 237 IF( sbc_blk_ice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard ice arrays' ) 238 #endif 192 239 ! 193 240 ! !** read bulk namelist … … 215 262 nblk = np_ECMWF ; ioptio = ioptio + 1 216 263 ENDIF 264 IF( ln_ANDREAS ) THEN 265 nblk = np_ANDREAS ; ioptio = ioptio + 1 266 ENDIF 217 267 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 218 268 … … 222 272 IF( ln_NCAR ) & 223 273 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' ) 274 IF( ln_ANDREAS ) & 275 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with ANDREAS algorithm' ) 224 276 IF( nn_fsbc /= 1 ) & 225 277 & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.') … … 254 306 ENDIF 255 307 ENDIF 308 309 #if defined key_si3 310 ioptio = 0 311 IF( ln_Cx_ice_cst ) THEN 312 nblk_ice = np_ice_cst ; ioptio = ioptio + 1 313 ENDIF 314 IF( ln_Cx_ice_LU12 ) THEN 315 nblk_ice = np_ice_lu12 ; ioptio = ioptio + 1 316 ENDIF 317 IF( ln_Cx_ice_LG15 ) THEN 318 nblk_ice = np_ice_lg15 ; ioptio = ioptio + 1 319 ENDIF 320 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one ice-atm bulk algorithm' ) 321 #endif 322 323 256 324 ! !* set the bulk structure 257 325 ! !- store namelist information in an array … … 310 378 ! 311 379 IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 ) & 312 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', &313 & ' This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' )380 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & 381 & ' This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' ) 314 382 ENDIF 315 383 END DO … … 321 389 !drag coefficient read from wave model definable only with mfs bulk formulae and core 322 390 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN 323 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and COREbulk formulae')391 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR bulk formulae') 324 392 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN 325 393 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') … … 341 409 ENDIF 342 410 ! 343 ! set transfer coefficients to default sea-ice values344 Cd_ice(:,:) = rCd_ice345 Ch_ice(:,:) = rCd_ice346 Ce_ice(:,:) = rCd_ice347 411 ! 348 412 IF(lwp) THEN !** Control print … … 350 414 WRITE(numout,*) !* namelist 351 415 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 352 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR416 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 353 417 WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0 354 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 355 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 418 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013) ln_COARE_3p6 = ', ln_COARE_3p6 419 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 420 WRITE(numout,*) ' "ANDREAS" algorithm (Andreas et al. 2015) ln_ANDREAS = ', ln_ANDREAS 356 421 WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt 357 422 WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu … … 359 424 WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac 360 425 WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))' 361 WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12362 WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15363 426 WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk 364 427 IF(ln_crt_fbk) THEN 365 WRITE(numout,*) ' Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta'366 WRITE(numout,*) ' Alpha rn_stau_a = ', rn_stau_a367 WRITE(numout,*) ' Beta rn_stau_b = ', rn_stau_b428 WRITE(numout,*) ' Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta' 429 WRITE(numout,*) ' Alpha rn_stau_a = ', rn_stau_a 430 WRITE(numout,*) ' Beta rn_stau_b = ', rn_stau_b 368 431 ENDIF 369 432 ! … … 374 437 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 375 438 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)' 439 CASE( np_ANDREAS ) ; WRITE(numout,*) ' ==>>> "ANDREAS" algorithm (Andreas et al. 2015)' 376 440 END SELECT 377 441 ! … … 386 450 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 387 451 END SELECT 452 ! 453 !#LB: 454 #if defined key_si3 455 IF( nn_ice > 0 ) THEN 456 WRITE(numout,*) 457 WRITE(numout,*) ' use constant ice-atm bulk transfer coeff. ln_Cx_ice_cst = ', ln_Cx_ice_cst 458 WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes, 2012 ln_Cx_ice_LU12 = ', ln_Cx_ice_LU12 459 WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes & Gryanik, 2015 ln_Cx_ice_LG15 = ', ln_Cx_ice_LG15 460 ENDIF 461 WRITE(numout,*) 462 SELECT CASE( nblk_ice ) !* Print the choice of bulk algorithm 463 CASE( np_ice_cst ) 464 WRITE(numout,*) ' ==>>> Constant bulk transfer coefficients over sea-ice:' 465 WRITE(numout,*) ' => Cd_ice, Ce_ice, Ch_ice =', REAL(rn_Cd_i,4), REAL(rn_Ce_i,4), REAL(rn_Ch_i,4) 466 IF( (rn_Cd_i<0._wp).OR.(rn_Cd_i>1.E-2_wp).OR.(rn_Ce_i<0._wp).OR.(rn_Ce_i>1.E-2_wp).OR.(rn_Ch_i<0._wp).OR.(rn_Ch_i>1.E-2_wp) ) & 467 & CALL ctl_stop( 'Be realistic in your pick of Cd_ice, Ce_ice & Ch_ice ! (0 < Cx < 1.E-2)') 468 CASE( np_ice_lu12 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes et al, 2012' 469 CASE( np_ice_lg15 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes & Gryanik, 2015' 470 END SELECT 471 #endif 472 !#LB. 388 473 ! 389 474 ENDIF … … 428 513 INTEGER, INTENT(in) :: kt ! ocean time step 429 514 !!---------------------------------------------------------------------- 430 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, z evp515 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp 431 516 REAL(wp) :: ztmp 432 517 !!---------------------------------------------------------------------- … … 465 550 ! ! compute the surface ocean fluxes using bulk formulea 466 551 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 552 553 ! Specific humidity of air at z=rn_zqt ! 554 SELECT CASE( nhumi ) 555 CASE( np_humi_sph ) 556 q_air_zt(:,:) = sf(jp_humi )%fnow(:,:,1) ! what we read in file is already a spec. humidity! 557 CASE( np_humi_dpt ) 558 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing q_air out of dew-point and P !' 559 q_air_zt(:,:) = q_sat( sf(jp_humi )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) 560 CASE( np_humi_rlh ) 561 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing q_air out of RH, t_air and slp !' !LBrm 562 q_air_zt(:,:) = q_air_rh( 0.01_wp*sf(jp_humi )%fnow(:,:,1), & 563 & sf(jp_tair )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) !#LB: 0.01 => RH is % percent in file 564 END SELECT 565 566 ! POTENTIAL temperature of air at z=rn_zqt 567 ! based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 568 ! (most reanalysis products provide absolute temp., not potential temp.) 569 IF( ln_tpot ) THEN 570 ! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...) 571 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing air POTENTIAL temperature out of ABSOLUTE temperature!' 572 theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt 573 ELSE 574 ! temperature read into file is already potential temperature 575 theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) 576 ENDIF 577 ! 467 578 CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1), & ! <<= in 468 & sf(jp_tair )%fnow(:,:,1), sf(jp_humi )%fnow(:,:,1),& ! <<= in579 & theta_air_zt(:,:), q_air_zt(:,:), & ! <<= in 469 580 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 470 581 & sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in 471 582 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 472 & tsk_m, zssq, zcd_du, zsen, z evp )! =>> out473 474 CALL blk_oce_2( sf(jp_tair )%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),& ! <<= in583 & tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out 584 585 CALL blk_oce_2( theta_air_zt(:,:), & ! <<= in 475 586 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1), & ! <<= in 476 587 & sf(jp_snow )%fnow(:,:,1), tsk_m, & ! <<= in 477 & zsen, z evp )! <=> in out588 & zsen, zlat, zevp ) ! <=> in out 478 589 ENDIF 479 590 ! … … 486 597 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 487 598 ENDIF 488 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 489 490 SELECT CASE( nhumi ) 491 CASE( np_humi_sph ) 492 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 493 CASE( np_humi_dpt ) 494 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 495 CASE( np_humi_rlh ) 496 qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file 497 END SELECT 599 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature instead ???? 600 !tatm_ice(:,:) = theta_air_zt(:,:) !#LB: THIS! ? 601 602 qatm_ice(:,:) = q_air_zt(:,:) !#LB: 498 603 499 604 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac … … 507 612 508 613 509 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, p humi, & ! inp614 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! inp 510 615 & pslp , pst , pu , pv, & ! inp 511 & puatm, pvatm, p qsr , pqlw ,& ! inp512 & ptsk , pssq , pcd_du, psen, p evp )! out616 & puatm, pvatm, pdqsr , pdqlw , & ! inp 617 & ptsk , pssq , pcd_du, psen, plat, pevp ) ! out 513 618 !!--------------------------------------------------------------------- 514 619 !! *** ROUTINE blk_oce_1 *** … … 523 628 !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) 524 629 !! - pcd_du : Cd x |dU| at T-points (m/s) 525 !! - psen : Ch x |dU| at T-points (m/s) 526 !! - pevp : Ce x |dU| at T-points (m/s) 630 !! - psen : sensible heat flux (W/m^2) 631 !! - plat : latent heat flux (W/m^2) 632 !! - pevp : evaporation (mm/s) #lolo 527 633 !!--------------------------------------------------------------------- 528 634 INTEGER , INTENT(in ) :: kt ! time step index 529 635 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s] 530 636 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s] 531 REAL(wp), INTENT(in ), DIMENSION(:,:) :: p humi! specific humidity at T-points [kg/kg]637 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqair ! specific humidity at T-points [kg/kg] 532 638 REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin] 533 639 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa] … … 537 643 REAL(wp), INTENT(in ), DIMENSION(:,:) :: puatm ! surface current seen by the atm at T-point (i-component) [m/s] 538 644 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pvatm ! surface current seen by the atm at T-point (j-component) [m/s] 539 REAL(wp), INTENT(in ), DIMENSION(:,:) :: p qsr !540 REAL(wp), INTENT(in ), DIMENSION(:,:) :: p qlw !645 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqsr ! downwelling solar (shortwave) radiation at surface [W/m^2] 646 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2] 541 647 REAL(wp), INTENT( out), DIMENSION(:,:) :: ptsk ! skin temp. (or SST if CS & WL not used) [Celsius] 542 648 REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] 543 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] 544 REAL(wp), INTENT( out), DIMENSION(:,:) :: psen ! Ch x |dU| at T-points [m/s] 545 REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp ! Ce x |dU| at T-points [m/s] 649 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du 650 REAL(wp), INTENT( out), DIMENSION(:,:) :: psen 651 REAL(wp), INTENT( out), DIMENSION(:,:) :: plat 652 REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp 546 653 ! 547 654 INTEGER :: ji, jj ! dummy loop indices … … 553 660 REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point 554 661 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 555 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K]556 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg]557 662 REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean 558 663 REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean 559 664 REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean 560 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux561 665 REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2 562 666 !!--------------------------------------------------------------------- … … 579 683 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 580 684 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 581 582 583 584 685 zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 686 zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 687 ! ... scalar wind at T-point (not masked) 688 wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) 585 689 END_2D 586 690 #else 587 691 ! ... scalar wind module at T-point (not masked) 588 692 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 589 693 wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 590 694 END_2D 591 695 #endif … … 597 701 zztmp = 1. - albo 598 702 IF( ln_dm2dc ) THEN 599 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)703 qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1) 600 704 ELSE 601 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)705 qsr(:,:) = zztmp * pdqsr(:,:) * tmask(:,:,1) 602 706 ENDIF 603 707 … … 616 720 ENDIF 617 721 618 ! specific humidity of air at "rn_zqt" m above the sea619 SELECT CASE( nhumi )620 CASE( np_humi_sph )621 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity!622 CASE( np_humi_dpt )623 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm624 zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) )625 CASE( np_humi_rlh )626 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm627 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file628 END SELECT629 !630 ! potential temperature of air at "rn_zqt" m above the sea631 IF( ln_abl ) THEN632 ztpot = ptair(:,:)633 ELSE634 ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate635 ! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2636 ! (since reanalysis products provide T at z, not theta !)637 !#LB: because AGRIF hates functions that return something else than a scalar, need to638 ! use scalar version of gamma_moist() ...639 IF( ln_tpot ) THEN640 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )641 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt642 END_2D643 ELSE644 ztpot = ptair(:,:)645 ENDIF646 ENDIF647 648 722 !! Time to call the user-selected bulk parameterization for 649 723 !! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more... … … 651 725 652 726 CASE( np_NCAR ) 653 CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, & 654 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 655 727 CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 728 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , & 729 & nb_iter=nn_iter_algo ) 730 ! 656 731 CASE( np_COARE_3p0 ) 657 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 658 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 659 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 660 732 CALL turb_coare3p0( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 733 & ln_skin_cs, ln_skin_wl, & 734 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 735 & nb_iter=nn_iter_algo, & 736 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 737 ! 661 738 CASE( np_COARE_3p6 ) 662 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 663 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 664 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 665 739 CALL turb_coare3p6( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 740 & ln_skin_cs, ln_skin_wl, & 741 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 742 & nb_iter=nn_iter_algo, & 743 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 744 ! 666 745 CASE( np_ECMWF ) 667 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 668 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 669 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 670 746 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 747 & ln_skin_cs, ln_skin_wl, & 748 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 749 & nb_iter=nn_iter_algo, & 750 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 751 ! 752 CASE( np_ANDREAS ) 753 CALL turb_andreas ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 754 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , & 755 & nb_iter=nn_iter_algo ) 756 ! 671 757 CASE DEFAULT 672 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formulaselected' )673 758 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' ) 759 ! 674 760 END SELECT 675 761 676 762 IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1)) 677 763 IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1)) 678 764 IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1)) 679 765 !! LB: mainly here for debugging purpose: 680 IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", ( ztpot-rt0) * tmask(:,:,1)) ! potential temperature at z=zt681 IF( iom_use('q_zt') ) CALL iom_put("q_zt", zqair * tmask(:,:,1)) ! specific humidity "682 IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (t _zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu766 IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ptair-rt0) * tmask(:,:,1)) ! potential temperature at z=zt 767 IF( iom_use('q_zt') ) CALL iom_put("q_zt", pqair * tmask(:,:,1)) ! specific humidity " 768 IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (theta_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu 683 769 IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity " 684 770 IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0 685 771 IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu 686 772 687 773 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 688 774 !! ptsk and pssq have been updated!!! … … 696 782 END IF 697 783 698 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbc blk_phy.F90784 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbc_phy.F90 699 785 ! ------------------------------------------------------------- 700 786 701 787 IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 702 788 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 703 704 705 706 707 708 rhoa(ji,jj) = rho_air( ptair(ji,jj), phumi(ji,jj), pslp(ji,jj) )789 zztmp = zU_zu(ji,jj) 790 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 791 pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 792 psen(ji,jj) = zztmp * zch_oce(ji,jj) 793 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 794 rhoa(ji,jj) = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) ) 709 795 END_2D 710 796 ELSE !== BLK formulation ==! turbulent fluxes computation 711 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t _zu(:,:), q_zu(:,:), &797 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), & 712 798 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 713 799 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 714 & taum(:,:), psen(:,:), zqla(:,:), &800 & taum(:,:), psen(:,:), plat(:,:), & 715 801 & pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac ) 716 802 717 zqla(:,:) = zqla(:,:) * tmask(:,:,1)718 803 psen(:,:) = psen(:,:) * tmask(:,:,1) 804 plat(:,:) = plat(:,:) * tmask(:,:,1) 719 805 taum(:,:) = taum(:,:) * tmask(:,:,1) 720 806 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 721 807 722 808 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 723 724 809 IF( wndm(ji,jj) > 0._wp ) THEN 810 zztmp = taum(ji,jj) / wndm(ji,jj) 725 811 #if defined key_cyclone 726 727 812 ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj) 813 ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj) 728 814 #else 729 730 731 #endif 732 733 734 ztau_j(ji,jj) = 0._wp735 815 ztau_i(ji,jj) = zztmp * pwndi(ji,jj) 816 ztau_j(ji,jj) = zztmp * pwndj(ji,jj) 817 #endif 818 ELSE 819 ztau_i(ji,jj) = 0._wp 820 ztau_j(ji,jj) = 0._wp 821 ENDIF 736 822 END_2D 737 823 … … 739 825 zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp ) ! set the max value of Stau corresponding to a wind of 3 m/s (<0) 740 826 DO_2D( 0, 1, 0, 1 ) ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop 741 742 743 744 827 zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) ! stau (<0) must be smaller than zstmax 828 ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) ) 829 ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) ) 830 taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) ) 745 831 END_2D 746 832 ENDIF … … 750 836 ! Note that coastal wind stress is not used in the code... so this extra care has no effect 751 837 DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T 752 753 754 755 838 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj ) ) & 839 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 840 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji ,jj+1) ) & 841 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 756 842 END_2D 843 757 844 758 845 IF( ln_crt_fbk ) THEN … … 762 849 ENDIF 763 850 764 CALL iom_put( "taum_oce", taum ) ! output wind stress module851 CALL iom_put( "taum_oce", taum*tmask(:,:,1) ) ! output wind stress module 765 852 766 853 IF(sn_cfctl%l_prtctl) THEN 767 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') 768 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & 769 & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 854 CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ') 855 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') 856 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & 857 & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 858 CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ') 770 859 ENDIF 771 860 ! 772 861 ENDIF !IF( ln_abl ) 773 862 774 863 ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius 775 864 776 865 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 777 866 CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius 778 867 CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference... 779 868 ENDIF 780 781 IF(sn_cfctl%l_prtctl) THEN782 CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' )783 CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' )784 CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' )785 ENDIF786 869 ! 787 870 END SUBROUTINE blk_oce_1 788 871 789 790 SUBROUTINE blk_oce_2( ptair, p qsr, pqlw, pprec,& ! <<= in791 & psnow, ptsk, psen, pevp ) ! <<= in872 873 SUBROUTINE blk_oce_2( ptair, pdqlw, pprec, psnow, & ! <<= in 874 & ptsk, psen, plat, pevp ) ! <<= in 792 875 !!--------------------------------------------------------------------- 793 876 !! *** ROUTINE blk_oce_2 *** … … 805 888 !! - emp : evaporation minus precipitation (kg/m2/s) 806 889 !!--------------------------------------------------------------------- 807 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair 808 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr 809 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw 890 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair ! potential temperature of air #LB: confirm! 891 REAL(wp), INTENT(in), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2] 810 892 REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec 811 893 REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow 812 894 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptsk ! SKIN surface temperature [Celsius] 813 895 REAL(wp), INTENT(in), DIMENSION(:,:) :: psen 896 REAL(wp), INTENT(in), DIMENSION(:,:) :: plat 814 897 REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp 815 898 ! 816 899 INTEGER :: ji, jj ! dummy loop indices 817 900 REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable 818 REAL(wp), DIMENSION(jpi,jpj) :: ztskk ! skin temp. in Kelvin 819 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes 820 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation 901 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! net long wave radiative heat flux 821 902 !!--------------------------------------------------------------------- 822 903 ! 823 904 ! local scalars ( place there for vector optimisation purposes) 824 905 825 826 ztskk(:,:) = ptsk(:,:) + rt0 ! => ptsk in Kelvin rather than Celsius827 828 906 ! ----------------------------------------------------------------------------- ! 829 907 ! III Net longwave radiative FLUX ! 830 908 ! ----------------------------------------------------------------------------- ! 831 832 !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 833 !! (ztskk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 834 zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*ztskk(:,:)*ztskk(:,:)*ztskk(:,:)*ztskk(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux 835 836 ! Latent flux over ocean 837 ! ----------------------- 838 839 ! use scalar version of L_vap() for AGRIF compatibility 840 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 841 zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj) ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 842 END_2D 843 844 IF(sn_cfctl%l_prtctl) THEN 845 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' ) 846 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 847 848 ENDIF 909 !! #LB: now moved after Turbulent fluxes because must use the skin temperature rather than bulk SST 910 !! (ptsk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 911 zqlw(:,:) = qlw_net( pdqlw(:,:), ptsk(:,:)+rt0 ) 849 912 850 913 ! ----------------------------------------------------------------------------- ! … … 855 918 & - pprec(:,:) * rn_pfac ) * tmask(:,:,1) 856 919 ! 857 qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar920 qns(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) & ! Downward Non Solar 858 921 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 859 922 & - pevp(:,:) * ptsk(:,:) * rcp & ! remove evap heat content at SST … … 865 928 ! 866 929 #if defined key_si3 867 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3)930 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) ! non solar without emp (only needed by SI3) 868 931 qsr_oce(:,:) = qsr(:,:) 869 932 #endif … … 873 936 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 874 937 CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean 875 CALL iom_put( "qla_oce" , zqla) ! output downward latent heat over the ocean938 CALL iom_put( "qla_oce" , plat ) ! output downward latent heat over the ocean 876 939 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 877 940 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] … … 880 943 ! 881 944 IF ( nn_ice == 0 ) THEN 882 CALL iom_put( "qemp_oce" , qns-zqlw-psen- zqla) ! output downward heat content of E-P over the ocean945 CALL iom_put( "qemp_oce" , qns-zqlw-psen-plat ) ! output downward heat content of E-P over the ocean 883 946 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 884 947 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean … … 888 951 IF(sn_cfctl%l_prtctl) THEN 889 952 CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ') 890 CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') 953 CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen : ' ) 954 CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat : ' ) 955 CALL prt_ctl(tab2d_1=qns , clinfo1=' blk_oce_2: qns : ' ) 891 956 CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ') 892 957 ENDIF … … 902 967 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 903 968 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 904 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag905 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag906 969 !!---------------------------------------------------------------------- 907 970 908 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, p humi, pslp , puice, pvice, ptsui, & ! inputs971 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, pqair, pslp , puice, pvice, ptsui, & ! inputs 909 972 & putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs 910 973 !!--------------------------------------------------------------------- … … 921 984 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s] 922 985 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s] 923 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: p humi! atmospheric wind at T-point [m/s]986 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric wind at T-point [m/s] 924 987 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s] 925 988 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! " … … 934 997 INTEGER :: ji, jj ! dummy loop indices 935 998 REAL(wp) :: zootm_su ! sea-ice surface mean temperature 936 REAL(wp) :: zztmp1, zztmp2 ! temporary arrays 937 REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau) 938 !!--------------------------------------------------------------------- 939 ! 940 999 REAL(wp) :: zztmp1, zztmp2 ! temporary scalars 1000 REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! temporary array 1001 !!--------------------------------------------------------------------- 1002 ! 1003 ! LB: ptsui is in K !!! 1004 ! 941 1005 ! ------------------------------------------------------------ ! 942 1006 ! Wind module relative to the moving ice ( U10m - U_ice ) ! … … 944 1008 ! C-grid ice dynamics : U & V-points (same as ocean) 945 1009 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 946 1010 wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 947 1011 END_2D 948 1012 ! 949 1013 ! Make ice-atm. drag dependent on ice concentration 950 IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations 951 CALL Cdn10_Lupkes2012( Cd_ice ) 952 Ch_ice(:,:) = Cd_ice(:,:) ! momentum and heat transfer coef. are considered identical 953 Ce_ice(:,:) = Cd_ice(:,:) 954 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 955 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 956 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 957 ENDIF 958 959 IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice) 960 IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice) 961 IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice) 962 963 ! local scalars ( place there for vector optimisation purposes) 964 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 1014 1015 1016 SELECT CASE( nblk_ice ) 1017 1018 CASE( np_ice_cst ) 1019 ! Constant bulk transfer coefficients over sea-ice: 1020 Cd_ice(:,:) = rn_Cd_i 1021 Ch_ice(:,:) = rn_Ch_i 1022 Ce_ice(:,:) = rn_Ce_i 1023 ! no height adjustment, keeping zt values: 1024 theta_zu_i(:,:) = ptair(:,:) 1025 q_zu_i(:,:) = pqair(:,:) 1026 1027 CASE( np_ice_lu12 ) 1028 ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 1029 CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, & 1030 & Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 1031 !! 1032 CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations 1033 ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 1034 CALL turb_ice_lg15( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, & 1035 & Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 1036 !! 1037 END SELECT 1038 1039 IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ai') ) & 1040 & ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice ! 1041 1042 IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp) 1043 IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp) 1044 IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp) 1045 965 1046 966 1047 IF( ln_blk ) THEN … … 969 1050 ! ---------------------------------------------------- ! 970 1051 ! supress moving ice in wind stress computation as we don't know how to do it properly... 971 DO_2D( 0, 1, 0, 1 ) ! at T point 972 putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndi(ji,jj) 973 pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndj(ji,jj) 1052 DO_2D( 0, 1, 0, 1 ) ! at T point 1053 zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj) 1054 putaui(ji,jj) = zztmp1 * pwndi(ji,jj) 1055 pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj) 974 1056 END_2D 1057 !#LB: saving the module of the ai wind-stress: NOT weighted by the ice concentration !!! 1058 IF(iom_use('taum_ai')) CALL iom_put( 'taum_ai', SQRT( putaui*putaui + pvtaui*pvtaui )*ztmp ) 975 1059 ! 976 1060 DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean). 977 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 978 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 979 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) 980 putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) ) 981 pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) ) 1061 !#LB: QUESTION?? so SI3 expects wind stress vector to be provided at U & V points? Not at T-points ? 1062 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 1063 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 1064 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) 1065 putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) ) 1066 pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) ) 982 1067 END_2D 983 1068 CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp ) … … 986 1071 & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) 987 1072 ELSE ! ln_abl 988 zztmp1 = 11637800.0_wp989 zztmp2 = -5897.8_wp990 1073 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 991 pcd_dui(ji,jj) = zcd_dui (ji,jj) 992 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 993 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 994 zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??) 995 pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 1074 pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj) 1075 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 1076 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 996 1077 END_2D 997 ENDIF 1078 !#LB: 1079 pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq ! 1080 !#LB. 1081 ENDIF !IF( ln_blk ) 998 1082 ! 999 1083 IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') … … 1002 1086 1003 1087 1004 SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, p humi, pslp, pqlw, pprec, psnow )1088 SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, pqair, pslp, pdqlw, pprec, psnow ) 1005 1089 !!--------------------------------------------------------------------- 1006 1090 !! *** ROUTINE blk_ice_2 *** … … 1018 1102 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 1019 1103 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 1020 REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair 1021 REAL(wp), DIMENSION(:,: ), INTENT(in) :: p humi1104 REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair ! potential temperature of air #LB: okay ??? 1105 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqair ! specific humidity of air 1022 1106 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp 1023 REAL(wp), DIMENSION(:,: ), INTENT(in) :: p qlw1107 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pdqlw 1024 1108 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec 1025 1109 REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow 1026 1110 !! 1027 1111 INTEGER :: ji, jj, jl ! dummy loop indices 1028 REAL(wp) :: zst 3! local variable1112 REAL(wp) :: zst, zst3, zsq ! local variable 1029 1113 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 1030 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 1031 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature 1114 REAL(wp) :: zztmp, zzblk, zztmp1, z1_rLsub ! - - 1032 1115 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice 1033 1116 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice … … 1035 1118 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 1036 1119 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3) 1037 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB1038 1120 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2 1039 1121 REAL(wp), DIMENSION(jpi,jpj) :: ztri 1040 1122 !!--------------------------------------------------------------------- 1041 1123 ! 1042 zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars 1043 zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 1044 ! 1045 SELECT CASE( nhumi ) 1046 CASE( np_humi_sph ) 1047 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 1048 CASE( np_humi_dpt ) 1049 zqair(:,:) = q_sat( phumi(:,:), pslp ) 1050 CASE( np_humi_rlh ) 1051 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 1052 END SELECT 1053 ! 1124 zcoef_dqlw = 4._wp * emiss_i * stefan ! local scalars 1125 ! 1126 1054 1127 zztmp = 1. / ( 1. - albo ) 1055 WHERE( ptsu(:,:,:) /= 0._wp ) 1056 z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 1057 ELSEWHERE 1058 z1_st(:,:,:) = 0._wp 1059 END WHERE 1128 dqla_ice(:,:,:) = 0._wp 1129 1060 1130 ! ! ========================== ! 1061 1131 DO jl = 1, jpl ! Loop over ice categories ! 1062 1132 ! ! ========================== ! 1063 DO jj = 1 , jpj 1064 DO ji = 1, jpi 1065 ! ----------------------------! 1066 ! I Radiative FLUXES ! 1067 ! ----------------------------! 1068 zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl) 1069 ! Short Wave (sw) 1070 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 1071 ! Long Wave (lw) 1072 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 1073 ! lw sensitivity 1074 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 1075 1076 ! ----------------------------! 1077 ! II Turbulent FLUXES ! 1078 ! ----------------------------! 1079 1080 ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 1081 ! Sensible Heat 1082 z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj)) 1083 ! Latent Heat 1084 zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 1085 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1086 & ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 1087 ! Latent heat sensitivity for ice (Dqla/Dt) 1088 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 1089 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1090 & z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2 1091 ELSE 1092 dqla_ice(ji,jj,jl) = 0._wp 1093 ENDIF 1094 1095 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 1096 z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 1097 1098 ! ----------------------------! 1099 ! III Total FLUXES ! 1100 ! ----------------------------! 1101 ! Downward Non Solar flux 1102 qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 1103 ! Total non solar heat flux sensitivity for ice 1104 dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) 1105 END DO 1106 ! 1107 END DO 1133 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1134 1135 zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K] 1136 zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present 1137 1138 ! ----------------------------! 1139 ! I Radiative FLUXES ! 1140 ! ----------------------------! 1141 ! Short Wave (sw) 1142 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 1143 1144 ! Long Wave (lw) 1145 zst3 = zst * zst * zst 1146 z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1) 1147 ! lw sensitivity 1148 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 1149 1150 ! ----------------------------! 1151 ! II Turbulent FLUXES ! 1152 ! ----------------------------! 1153 1154 ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 1155 1156 ! Common term in bulk F. equations... 1157 zzblk = rhoa(ji,jj) * wndm_ice(ji,jj) 1158 1159 ! Sensible Heat 1160 zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj) 1161 z_qsb (ji,jj,jl) = zztmp1 * (zst - theta_zu_i(ji,jj)) 1162 z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice) 1163 1164 ! Latent Heat 1165 zztmp1 = zzblk * rLsub * Ce_ice(ji,jj) 1166 qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ??? 1167 IF(qla_ice(ji,jj,jl)>0._wp) dqla_ice(ji,jj,jl) = zztmp1*dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT) 1168 ! !#LB: dq_sat_dt_ice() in "sbc_phy.F90" 1169 !#LB: without this unjustified "condensation sensure": 1170 !qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj)) 1171 !dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT) 1172 1173 1174 1175 ! ----------------------------! 1176 ! III Total FLUXES ! 1177 ! ----------------------------! 1178 1179 ! Downward Non Solar flux 1180 qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 1181 1182 ! Total non solar heat flux sensitivity for ice 1183 dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ???? 1184 1185 END_2D 1108 1186 ! 1109 1187 END DO … … 1157 1235 ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm 1158 1236 DO jl = 1, jpl 1159 WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1237 WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1160 1238 qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 1161 1239 ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm 1162 1240 qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:) 1163 1241 ELSEWHERE ! zero when hs>0 1164 qtr_ice_top(:,:,jl) = 0._wp 1242 qtr_ice_top(:,:,jl) = 0._wp 1165 1243 END WHERE 1166 1244 ENDDO … … 1201 1279 CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ') 1202 1280 ENDIF 1203 ! 1281 1282 !#LB: 1283 ! air-ice heat flux components that are not written from ice_stp()@icestp.F90: 1284 IF( iom_use('qla_ice') ) CALL iom_put( 'qla_ice', SUM( - qla_ice * a_i_b, dim=3 ) ) !#LB: sign consistent with what's done for ocean 1285 IF( iom_use('qsb_ice') ) CALL iom_put( 'qsb_ice', SUM( - z_qsb * a_i_b, dim=3 ) ) !#LB: ==> negative => loss of heat for sea-ice 1286 IF( iom_use('qlw_ice') ) CALL iom_put( 'qlw_ice', SUM( z_qlw * a_i_b, dim=3 ) ) 1287 !#LB. 1288 1204 1289 END SUBROUTINE blk_ice_2 1205 1290 … … 1254 1339 DO jl = 1, jpl 1255 1340 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1256 1257 1341 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness 1342 IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 1258 1343 END_2D 1259 1344 END DO … … 1269 1354 DO jl = 1, jpl 1270 1355 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1356 ! 1357 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1358 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 1359 ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature 1360 ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature 1361 zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 1362 ! 1363 DO iter = 1, nit ! --- Iterative loop 1364 zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards) 1365 zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget 1366 ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update 1367 END DO 1368 ! 1369 ptsu (ji,jj,jl) = MIN( rt0, ztsu ) 1370 qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 1371 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1372 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) ) & 1373 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1374 1375 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1376 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) 1292 1377 1293 1378 END_2D … … 1296 1381 ! 1297 1382 END SUBROUTINE blk_ice_qcn 1298 1299 1300 SUBROUTINE Cdn10_Lupkes2012( pcd )1301 !!----------------------------------------------------------------------1302 !! *** ROUTINE Cdn10_Lupkes2012 ***1303 !!1304 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m1305 !! to make it dependent on edges at leads, melt ponds and flows.1306 !! After some approximations, this can be resumed to a dependency1307 !! on ice concentration.1308 !!1309 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50)1310 !! with the highest level of approximation: level4, eq.(59)1311 !! The generic drag over a cell partly covered by ice can be re-written as follows:1312 !!1313 !! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu1314 !!1315 !! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59)1316 !! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)1317 !! A is the concentration of ice minus melt ponds (if any)1318 !!1319 !! This new drag has a parabolic shape (as a function of A) starting at1320 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.51321 !! and going down to Cdi(say 1.4e-3) for A=11322 !!1323 !! It is theoretically applicable to all ice conditions (not only MIZ)1324 !! => see Lupkes et al (2013)1325 !!1326 !! ** References : Lupkes et al. JGR 2012 (theory)1327 !! Lupkes et al. GRL 2013 (application to GCM)1328 !!1329 !!----------------------------------------------------------------------1330 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd1331 REAL(wp), PARAMETER :: zCe = 2.23e-03_wp1332 REAL(wp), PARAMETER :: znu = 1._wp1333 REAL(wp), PARAMETER :: zmu = 1._wp1334 REAL(wp), PARAMETER :: zbeta = 1._wp1335 REAL(wp) :: zcoef1336 !!----------------------------------------------------------------------1337 zcoef = znu + 1._wp / ( 10._wp * zbeta )1338 1339 ! generic drag over a cell partly covered by ice1340 !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag1341 !! & Cd_ice * at_i_b(:,:) + & ! pure ice drag1342 !! & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu ! change due to sea-ice morphology1343 1344 ! ice-atm drag1345 pcd(:,:) = rCd_ice + & ! pure ice drag1346 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology1347 1348 END SUBROUTINE Cdn10_Lupkes20121349 1350 1351 SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch )1352 !!----------------------------------------------------------------------1353 !! *** ROUTINE Cdn10_Lupkes2015 ***1354 !!1355 !! ** pUrpose : Alternative turbulent transfert coefficients formulation1356 !! between sea-ice and atmosphere with distinct momentum1357 !! and heat coefficients depending on sea-ice concentration1358 !! and atmospheric stability (no meltponds effect for now).1359 !!1360 !! ** Method : The parameterization is adapted from Lupkes et al. (2015)1361 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,1362 !! it considers specific skin and form drags (Andreas et al. 2010)1363 !! to compute neutral transfert coefficients for both heat and1364 !! momemtum fluxes. Atmospheric stability effect on transfert1365 !! coefficient is also taken into account following Louis (1979).1366 !!1367 !! ** References : Lupkes et al. JGR 2015 (theory)1368 !! Lupkes et al. ECHAM6 documentation 2015 (implementation)1369 !!1370 !!----------------------------------------------------------------------1371 !1372 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K]1373 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa]1374 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient1375 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient1376 REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat1377 !1378 ! ECHAM6 constants1379 REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m]1380 REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m]1381 REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m]1382 REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 411383 REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 411384 REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 131385 REAL(wp), PARAMETER :: zc2 = zc * zc1386 REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 141387 REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 301388 REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 511389 REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 561390 REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 261391 REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 261392 REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma1393 REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp1394 !1395 INTEGER :: ji, jj ! dummy loop indices1396 REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu1397 REAL(wp) :: zrib_o, zrib_i1398 REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice1399 REAL(wp) :: zChn_skin_ice, zChn_form_ice1400 REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw1401 REAL(wp) :: zCdn_form_tmp1402 !!----------------------------------------------------------------------1403 1404 ! Momentum Neutral Transfert Coefficients (should be a constant)1405 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 401406 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 71407 zCdn_ice = zCdn_skin_ice ! Eq. 71408 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32)1409 1410 ! Heat Neutral Transfert Coefficients1411 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. 521412 1413 ! Atmospheric and Surface Variables1414 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin1415 zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg]1416 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg]1417 !1418 DO_2D( 0, 0, 0, 0 )1419 ! Virtual potential temperature [K]1420 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean1421 zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice1422 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu1423 1424 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)1425 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean1426 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice1427 1428 ! Momentum and Heat Neutral Transfert Coefficients1429 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 401430 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 531431 1432 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?)1433 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water1434 z0i = z0_skin_ice ! over ice1435 IF( zrib_o <= 0._wp ) THEN1436 zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 101437 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 261438 & )**zgamma )**z1_gamma1439 ELSE1440 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 121441 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 281442 ENDIF1443 1444 IF( zrib_i <= 0._wp ) THEN1445 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 91446 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 251447 ELSE1448 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 111449 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 271450 ENDIF1451 1452 ! Momentum Transfert Coefficients (Eq. 38)1453 pcd(ji,jj) = zCdn_skin_ice * zfmi + &1454 & 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) )1455 1456 ! Heat Transfert Coefficients (Eq. 49)1457 pch(ji,jj) = zChn_skin_ice * zfhi + &1458 & zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) )1459 !1460 END_2D1461 CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1.0_wp, pch, 'T', 1.0_wp )1462 !1463 END SUBROUTINE Cdn10_Lupkes20151464 1383 1465 1384 #endif
Note: See TracChangeset
for help on using the changeset viewer.
