Changeset 14072 for NEMO/trunk/src/OCE/SBC/sbcblk.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/SBC/sbcblk.F90 (modified) (51 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/SBC/sbcblk.F90
r14007 r14072 19 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 20 20 !! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE) 21 !! 4.2 ! 2020-12 (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements 21 22 !!---------------------------------------------------------------------- 22 23 … … 30 31 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 31 32 !! 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 33 !!---------------------------------------------------------------------- 35 34 USE oce ! ocean dynamics and tracers … … 41 40 USE sbcdcy ! surface boundary condition: diurnal cycle 42 41 USE sbcwave , ONLY : cdn_wave ! wave module 43 USE sbc_ice ! Surface boundary condition: ice fields44 42 USE lib_fortran ! to use key_nosignedzero 43 ! 45 44 #if defined key_si3 45 USE sbc_ice ! Surface boundary condition: ice fields #LB? ok to be in 'key_si3' ??? 46 46 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 47 USE icevar ! for CALL ice_var_snwblow 48 #endif 49 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 48 USE sbcblk_algo_ice_an05 49 USE sbcblk_algo_ice_lu12 50 USE sbcblk_algo_ice_lg15 51 #endif 52 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - (formerly known as CORE, Large & Yeager, 2009) 50 53 USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 51 54 USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 52 55 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 45r1) 56 USE sbcblk_algo_andreas ! => turb_andreas : Andreas et al. 2015 53 57 ! 54 58 USE iom ! I/O manager library … … 58 62 USE prtctl ! Print control 59 63 60 USE sbcblk_phy ! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc... 61 64 USE sbc_phy ! Catalog of functions for physical/meteorological parameters in the marine boundary layer 62 65 63 66 IMPLICIT NONE … … 100 103 LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013) 101 104 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1) 105 LOGICAL :: ln_ANDREAS ! "ANDREAS" algorithm (Andreas et al. 2015) 102 106 ! 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) 107 !#LB: 108 LOGICAL :: ln_Cx_ice_cst ! use constant air-ice bulk transfer coefficients (value given in namelist's rn_Cd_i, rn_Ce_i & rn_Ch_i) 109 REAL(wp) :: rn_Cd_i, rn_Ce_i, rn_Ch_i ! values for " " 110 LOGICAL :: ln_Cx_ice_AN05 ! air-ice bulk transfer coefficients based on Andreas et al., 2005 111 LOGICAL :: ln_Cx_ice_LU12 ! air-ice bulk transfer coefficients based on Lupkes et al., 2012 112 LOGICAL :: ln_Cx_ice_LG15 ! air-ice bulk transfer coefficients based on Lupkes & Gryanik, 2015 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 … … 122 137 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 123 138 LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB 124 LOGICAL :: ln_t pot !!GS: flag to compute or not potential temperature139 LOGICAL :: ln_tair_pot ! temperature read in files ("sn_tair") is already potential temperature (not absolute) 125 140 ! 126 141 INTEGER :: nhumi ! choice of the bulk algorithm … … 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_an05 = 2 ! Andreas et al., 2005 162 INTEGER, PARAMETER :: np_ice_lu12 = 3 ! Lupkes el al., 2012 163 INTEGER, PARAMETER :: np_ice_lg15 = 4 ! Lupkes & Gryanik, 2015 164 #endif 165 !LB. 166 167 138 168 139 169 !! * Substitutions … … 150 180 !! *** ROUTINE sbc_blk_alloc *** 151 181 !!------------------------------------------------------------------- 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 ! 182 ALLOCATE( theta_zu(jpi,jpj), q_zu(jpi,jpj), STAT=sbc_blk_alloc ) 156 183 CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) 157 184 IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' ) 158 185 END FUNCTION sbc_blk_alloc 186 187 #if defined key_si3 188 INTEGER FUNCTION sbc_blk_ice_alloc() 189 !!------------------------------------------------------------------- 190 !! *** ROUTINE sbc_blk_ice_alloc *** 191 !!------------------------------------------------------------------- 192 ALLOCATE( Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), & 193 & theta_zu_i(jpi,jpj), q_zu_i(jpi,jpj), STAT=sbc_blk_ice_alloc ) 194 CALL mpp_sum ( 'sbcblk', sbc_blk_ice_alloc ) 195 IF( sbc_blk_ice_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_ice_alloc: failed to allocate arrays' ) 196 END FUNCTION sbc_blk_ice_alloc 197 #endif 159 198 160 199 … … 178 217 TYPE(FLD_N) :: sn_cc, sn_hpgi, sn_hpgj ! " " 179 218 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, & 219 NAMELIST/namsbc_blk/ ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, ln_ANDREAS, & ! bulk algorithm 220 & rn_zqt, rn_zu, nn_iter_algo, ln_skin_cs, ln_skin_wl, & 221 & rn_pfac, rn_efac, & 186 222 & 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 223 & ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tair_pot, & 224 & ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i, & 225 & ln_Cx_ice_AN05, ln_Cx_ice_LU12, ln_Cx_ice_LG15, & 226 & cn_dir, & 227 & sn_wndi, sn_wndj, sn_qsr, sn_qlw , & ! input fields 228 & sn_tair, sn_humi, sn_prec, sn_snow, sn_slp, & 229 & sn_uoatm, sn_voatm, sn_cc, sn_hpgi, sn_hpgj 230 231 ! cool-skin / warm-layer !LB 188 232 !!--------------------------------------------------------------------- 189 233 ! 190 234 ! ! allocate sbc_blk_core array 191 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 235 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 236 ! 237 #if defined key_si3 238 IF( sbc_blk_ice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard ice arrays' ) 239 #endif 192 240 ! 193 241 ! !** read bulk namelist … … 215 263 nblk = np_ECMWF ; ioptio = ioptio + 1 216 264 ENDIF 265 IF( ln_ANDREAS ) THEN 266 nblk = np_ANDREAS ; ioptio = ioptio + 1 267 ENDIF 217 268 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 218 269 … … 222 273 IF( ln_NCAR ) & 223 274 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' ) 275 IF( ln_ANDREAS ) & 276 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with ANDREAS algorithm' ) 224 277 IF( nn_fsbc /= 1 ) & 225 278 & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.') … … 254 307 ENDIF 255 308 ENDIF 309 310 #if defined key_si3 311 ioptio = 0 312 IF( ln_Cx_ice_cst ) THEN 313 nblk_ice = np_ice_cst ; ioptio = ioptio + 1 314 ENDIF 315 IF( ln_Cx_ice_AN05 ) THEN 316 nblk_ice = np_ice_an05 ; ioptio = ioptio + 1 317 ENDIF 318 IF( ln_Cx_ice_LU12 ) THEN 319 nblk_ice = np_ice_lu12 ; ioptio = ioptio + 1 320 ENDIF 321 IF( ln_Cx_ice_LG15 ) THEN 322 nblk_ice = np_ice_lg15 ; ioptio = ioptio + 1 323 ENDIF 324 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one ice-atm bulk algorithm' ) 325 #endif 326 327 256 328 ! !* set the bulk structure 257 329 ! !- store namelist information in an array … … 322 394 ENDIF 323 395 ! 324 ! set transfer coefficients to default sea-ice values325 Cd_ice(:,:) = rCd_ice326 Ch_ice(:,:) = rCd_ice327 Ce_ice(:,:) = rCd_ice328 396 ! 329 397 IF(lwp) THEN !** Control print … … 331 399 WRITE(numout,*) !* namelist 332 400 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 333 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR401 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 334 402 WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0 335 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 336 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 403 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013) ln_COARE_3p6 = ', ln_COARE_3p6 404 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 405 WRITE(numout,*) ' "ANDREAS" algorithm (Andreas et al. 2015) ln_ANDREAS = ', ln_ANDREAS 337 406 WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt 338 407 WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu … … 340 409 WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac 341 410 WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))' 342 WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12343 WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15344 411 WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk 345 412 IF(ln_crt_fbk) THEN … … 355 422 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 356 423 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)' 424 CASE( np_ANDREAS ) ; WRITE(numout,*) ' ==>>> "ANDREAS" algorithm (Andreas et al. 2015)' 357 425 END SELECT 358 426 ! … … 367 435 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 368 436 END SELECT 437 ! 438 !#LB: 439 #if defined key_si3 440 IF( nn_ice > 0 ) THEN 441 WRITE(numout,*) 442 WRITE(numout,*) ' use constant ice-atm bulk transfer coeff. ln_Cx_ice_cst = ', ln_Cx_ice_cst 443 WRITE(numout,*) ' use ice-atm bulk coeff. from Andreas et al., 2005 ln_Cx_ice_AN05 = ', ln_Cx_ice_AN05 444 WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes et al., 2012 ln_Cx_ice_LU12 = ', ln_Cx_ice_LU12 445 WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes & Gryanik, 2015 ln_Cx_ice_LG15 = ', ln_Cx_ice_LG15 446 ENDIF 447 WRITE(numout,*) 448 SELECT CASE( nblk_ice ) !* Print the choice of bulk algorithm 449 CASE( np_ice_cst ) 450 WRITE(numout,*) ' ==>>> Constant bulk transfer coefficients over sea-ice:' 451 WRITE(numout,*) ' => Cd_ice, Ce_ice, Ch_ice =', REAL(rn_Cd_i,4), REAL(rn_Ce_i,4), REAL(rn_Ch_i,4) 452 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) ) & 453 & CALL ctl_stop( 'Be realistic in your pick of Cd_ice, Ce_ice & Ch_ice ! (0 < Cx < 1.E-2)') 454 CASE( np_ice_an05 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Andreas et al, 2005' 455 CASE( np_ice_lu12 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes et al, 2012' 456 CASE( np_ice_lg15 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes & Gryanik, 2015' 457 END SELECT 458 #endif 459 !#LB. 369 460 ! 370 461 ENDIF … … 409 500 INTEGER, INTENT(in) :: kt ! ocean time step 410 501 !!---------------------------------------------------------------------- 411 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, z evp502 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp 412 503 REAL(wp) :: ztmp 413 504 !!---------------------------------------------------------------------- … … 446 537 ! ! compute the surface ocean fluxes using bulk formulea 447 538 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 539 540 ! Specific humidity of air at z=rn_zqt ! 541 SELECT CASE( nhumi ) 542 CASE( np_humi_sph ) 543 q_air_zt(:,:) = sf(jp_humi )%fnow(:,:,1) ! what we read in file is already a spec. humidity! 544 CASE( np_humi_dpt ) 545 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing q_air out of dew-point and P !' 546 q_air_zt(:,:) = q_sat( sf(jp_humi )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) 547 CASE( np_humi_rlh ) 548 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing q_air out of RH, t_air and slp !' !LBrm 549 q_air_zt(:,:) = q_air_rh( 0.01_wp*sf(jp_humi )%fnow(:,:,1), & 550 & sf(jp_tair )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) !#LB: 0.01 => RH is % percent in file 551 END SELECT 552 553 ! POTENTIAL temperature of air at z=rn_zqt 554 ! based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 555 ! (most reanalysis products provide absolute temp., not potential temp.) 556 IF( ln_tair_pot ) THEN 557 ! temperature read into file is already potential temperature, do nothing... 558 theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) 559 ELSE 560 ! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...) 561 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!' 562 theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt 563 ENDIF 564 ! 448 565 CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1), & ! <<= in 449 & sf(jp_tair )%fnow(:,:,1), sf(jp_humi )%fnow(:,:,1),& ! <<= in566 & theta_air_zt(:,:), q_air_zt(:,:), & ! <<= in 450 567 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 451 568 & sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in 452 569 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 453 & tsk_m, zssq, zcd_du, zsen, z evp )! =>> out454 455 CALL blk_oce_2( sf(jp_tair )%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),& ! <<= in570 & tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out 571 572 CALL blk_oce_2( theta_air_zt(:,:), & ! <<= in 456 573 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1), & ! <<= in 457 574 & sf(jp_snow )%fnow(:,:,1), tsk_m, & ! <<= in 458 & zsen, z evp )! <=> in out575 & zsen, zlat, zevp ) ! <=> in out 459 576 ENDIF 460 577 ! … … 467 584 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 468 585 ENDIF 469 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 470 471 SELECT CASE( nhumi ) 472 CASE( np_humi_sph ) 473 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 474 CASE( np_humi_dpt ) 475 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 476 CASE( np_humi_rlh ) 477 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 478 END SELECT 586 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature instead ???? 587 !tatm_ice(:,:) = theta_air_zt(:,:) !#LB: THIS! ? 588 589 qatm_ice(:,:) = q_air_zt(:,:) !#LB: 479 590 480 591 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac … … 488 599 489 600 490 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, p humi, & ! inp601 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! inp 491 602 & pslp , pst , pu , pv, & ! inp 492 & puatm, pvatm, p qsr , pqlw ,& ! inp493 & ptsk , pssq , pcd_du, psen, p evp )! out603 & puatm, pvatm, pdqsr , pdqlw , & ! inp 604 & ptsk , pssq , pcd_du, psen, plat, pevp ) ! out 494 605 !!--------------------------------------------------------------------- 495 606 !! *** ROUTINE blk_oce_1 *** … … 504 615 !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) 505 616 !! - pcd_du : Cd x |dU| at T-points (m/s) 506 !! - psen : Ch x |dU| at T-points (m/s) 507 !! - pevp : Ce x |dU| at T-points (m/s) 617 !! - psen : sensible heat flux (W/m^2) 618 !! - plat : latent heat flux (W/m^2) 619 !! - pevp : evaporation (mm/s) #lolo 508 620 !!--------------------------------------------------------------------- 509 621 INTEGER , INTENT(in ) :: kt ! time step index 510 622 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s] 511 623 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s] 512 REAL(wp), INTENT(in ), DIMENSION(:,:) :: p humi! specific humidity at T-points [kg/kg]624 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqair ! specific humidity at T-points [kg/kg] 513 625 REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin] 514 626 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa] … … 518 630 REAL(wp), INTENT(in ), DIMENSION(:,:) :: puatm ! surface current seen by the atm at T-point (i-component) [m/s] 519 631 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pvatm ! surface current seen by the atm at T-point (j-component) [m/s] 520 REAL(wp), INTENT(in ), DIMENSION(:,:) :: p qsr !521 REAL(wp), INTENT(in ), DIMENSION(:,:) :: p qlw !632 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqsr ! downwelling solar (shortwave) radiation at surface [W/m^2] 633 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2] 522 634 REAL(wp), INTENT( out), DIMENSION(:,:) :: ptsk ! skin temp. (or SST if CS & WL not used) [Celsius] 523 635 REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] 524 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] 525 REAL(wp), INTENT( out), DIMENSION(:,:) :: psen ! Ch x |dU| at T-points [m/s] 526 REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp ! Ce x |dU| at T-points [m/s] 636 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du 637 REAL(wp), INTENT( out), DIMENSION(:,:) :: psen 638 REAL(wp), INTENT( out), DIMENSION(:,:) :: plat 639 REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp 527 640 ! 528 641 INTEGER :: ji, jj ! dummy loop indices … … 534 647 REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point 535 648 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 536 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K]537 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg]538 649 REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean 539 650 REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean 540 651 REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean 541 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux542 652 REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2 543 653 !!--------------------------------------------------------------------- … … 578 688 zztmp = 1. - albo 579 689 IF( ln_dm2dc ) THEN 580 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)690 qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1) 581 691 ELSE 582 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)692 qsr(:,:) = zztmp * pdqsr(:,:) * tmask(:,:,1) 583 693 ENDIF 584 694 … … 597 707 ENDIF 598 708 599 ! specific humidity of air at "rn_zqt" m above the sea600 SELECT CASE( nhumi )601 CASE( np_humi_sph )602 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity!603 CASE( np_humi_dpt )604 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm605 zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) )606 CASE( np_humi_rlh )607 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm608 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file609 END SELECT610 !611 ! potential temperature of air at "rn_zqt" m above the sea612 IF( ln_abl ) THEN613 ztpot = ptair(:,:)614 ELSE615 ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate616 ! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2617 ! (since reanalysis products provide T at z, not theta !)618 !#LB: because AGRIF hates functions that return something else than a scalar, need to619 ! use scalar version of gamma_moist() ...620 IF( ln_tpot ) THEN621 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )622 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt623 END_2D624 ELSE625 ztpot = ptair(:,:)626 ENDIF627 ENDIF628 629 709 !! Time to call the user-selected bulk parameterization for 630 710 !! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more... … … 632 712 633 713 CASE( np_NCAR ) 634 CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, & 635 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 636 714 CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 715 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , & 716 & nb_iter=nn_iter_algo ) 717 ! 637 718 CASE( np_COARE_3p0 ) 638 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 639 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 640 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 641 719 CALL turb_coare3p0( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 720 & ln_skin_cs, ln_skin_wl, & 721 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 722 & nb_iter=nn_iter_algo, & 723 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 724 ! 642 725 CASE( np_COARE_3p6 ) 643 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 644 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 645 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 646 726 CALL turb_coare3p6( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 727 & ln_skin_cs, ln_skin_wl, & 728 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 729 & nb_iter=nn_iter_algo, & 730 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 731 ! 647 732 CASE( np_ECMWF ) 648 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 649 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 650 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 651 733 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 734 & ln_skin_cs, ln_skin_wl, & 735 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 736 & nb_iter=nn_iter_algo, & 737 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 738 ! 739 CASE( np_ANDREAS ) 740 CALL turb_andreas ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 741 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , & 742 & nb_iter=nn_iter_algo ) 743 ! 652 744 CASE DEFAULT 653 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formulaselected' )654 745 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' ) 746 ! 655 747 END SELECT 656 748 657 749 IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1)) 658 750 IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1)) 659 751 IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1)) 660 752 !! LB: mainly here for debugging purpose: 661 IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", ( ztpot-rt0) * tmask(:,:,1)) ! potential temperature at z=zt662 IF( iom_use('q_zt') ) CALL iom_put("q_zt", zqair * tmask(:,:,1)) ! specific humidity "663 IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (t _zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu753 IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ptair-rt0) * tmask(:,:,1)) ! potential temperature at z=zt 754 IF( iom_use('q_zt') ) CALL iom_put("q_zt", pqair * tmask(:,:,1)) ! specific humidity " 755 IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (theta_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu 664 756 IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity " 665 757 IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0 666 758 IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu 667 759 668 760 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 669 761 !! ptsk and pssq have been updated!!! … … 677 769 END IF 678 770 679 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbc blk_phy.F90771 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbc_phy.F90 680 772 ! ------------------------------------------------------------- 681 773 … … 687 779 psen(ji,jj) = zztmp * zch_oce(ji,jj) 688 780 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 689 rhoa(ji,jj) = rho_air( ptair(ji,jj), p humi(ji,jj), pslp(ji,jj) )781 rhoa(ji,jj) = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) ) 690 782 END_2D 691 783 ELSE !== BLK formulation ==! turbulent fluxes computation 692 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t _zu(:,:), q_zu(:,:), &784 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), & 693 785 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 694 786 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 695 & taum(:,:), psen(:,:), zqla(:,:), &787 & taum(:,:), psen(:,:), plat(:,:), & 696 788 & pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac ) 697 789 698 zqla(:,:) = zqla(:,:) * tmask(:,:,1)699 790 psen(:,:) = psen(:,:) * tmask(:,:,1) 791 plat(:,:) = plat(:,:) * tmask(:,:,1) 700 792 taum(:,:) = taum(:,:) * tmask(:,:,1) 701 793 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 702 794 703 795 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 704 705 796 IF( wndm(ji,jj) > 0._wp ) THEN 797 zztmp = taum(ji,jj) / wndm(ji,jj) 706 798 #if defined key_cyclone 707 708 799 ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj) 800 ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj) 709 801 #else 710 711 712 #endif 713 714 715 ztau_j(ji,jj) = 0._wp716 802 ztau_i(ji,jj) = zztmp * pwndi(ji,jj) 803 ztau_j(ji,jj) = zztmp * pwndj(ji,jj) 804 #endif 805 ELSE 806 ztau_i(ji,jj) = 0._wp 807 ztau_j(ji,jj) = 0._wp 808 ENDIF 717 809 END_2D 718 810 … … 743 835 ENDIF 744 836 745 CALL iom_put( "taum_oce", taum ) ! output wind stress module 837 ! Saving open-ocean wind-stress (module and components) on T-points: 838 CALL iom_put( "taum_oce", taum(:,:)*tmask(:,:,1) ) ! output wind stress module 839 !#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau" (U-grid) and vtau" (V-grid) does the job in: [DYN/dynatf.F90]) 840 CALL iom_put( "utau_oce", ztau_i(:,:)*tmask(:,:,1) ) ! utau at T-points! 841 CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) ) ! vtau at T-points! 746 842 747 843 IF(sn_cfctl%l_prtctl) THEN 748 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') 749 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & 750 & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 844 CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ') 845 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') 846 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & 847 & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 848 CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ') 751 849 ENDIF 752 850 ! 753 851 ENDIF !IF( ln_abl ) 754 852 755 853 ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius 756 854 757 855 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 758 856 CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius 759 857 CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference... 760 858 ENDIF 761 762 IF(sn_cfctl%l_prtctl) THEN763 CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' )764 CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' )765 CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' )766 ENDIF767 859 ! 768 860 END SUBROUTINE blk_oce_1 769 861 770 862 771 SUBROUTINE blk_oce_2( ptair, p qsr, pqlw, pprec,& ! <<= in772 & psnow, ptsk, psen, pevp ) ! <<= in863 SUBROUTINE blk_oce_2( ptair, pdqlw, pprec, psnow, & ! <<= in 864 & ptsk, psen, plat, pevp ) ! <<= in 773 865 !!--------------------------------------------------------------------- 774 866 !! *** ROUTINE blk_oce_2 *** … … 786 878 !! - emp : evaporation minus precipitation (kg/m2/s) 787 879 !!--------------------------------------------------------------------- 788 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair 789 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr 790 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw 880 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair ! potential temperature of air #LB: confirm! 881 REAL(wp), INTENT(in), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2] 791 882 REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec 792 883 REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow 793 884 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptsk ! SKIN surface temperature [Celsius] 794 885 REAL(wp), INTENT(in), DIMENSION(:,:) :: psen 886 REAL(wp), INTENT(in), DIMENSION(:,:) :: plat 795 887 REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp 796 888 ! 797 889 INTEGER :: ji, jj ! dummy loop indices 798 890 REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable 799 REAL(wp), DIMENSION(jpi,jpj) :: ztskk ! skin temp. in Kelvin 800 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes 801 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation 891 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! net long wave radiative heat flux 802 892 !!--------------------------------------------------------------------- 803 893 ! 804 894 ! local scalars ( place there for vector optimisation purposes) 805 895 806 807 ztskk(:,:) = ptsk(:,:) + rt0 ! => ptsk in Kelvin rather than Celsius808 809 896 ! ----------------------------------------------------------------------------- ! 810 897 ! III Net longwave radiative FLUX ! 811 898 ! ----------------------------------------------------------------------------- ! 812 813 !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 814 !! (ztskk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 815 zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*ztskk(:,:)*ztskk(:,:)*ztskk(:,:)*ztskk(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux 816 817 ! Latent flux over ocean 818 ! ----------------------- 819 820 ! use scalar version of L_vap() for AGRIF compatibility 821 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 822 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 823 END_2D 824 825 IF(sn_cfctl%l_prtctl) THEN 826 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' ) 827 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 828 829 ENDIF 899 !! #LB: now moved after Turbulent fluxes because must use the skin temperature rather than bulk SST 900 !! (ptsk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 901 zqlw(:,:) = qlw_net( pdqlw(:,:), ptsk(:,:)+rt0 ) 830 902 831 903 ! ----------------------------------------------------------------------------- ! … … 836 908 & - pprec(:,:) * rn_pfac ) * tmask(:,:,1) 837 909 ! 838 qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar910 qns(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) & ! Downward Non Solar 839 911 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 840 912 & - pevp(:,:) * ptsk(:,:) * rcp & ! remove evap heat content at SST … … 846 918 ! 847 919 #if defined key_si3 848 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3)920 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) ! non solar without emp (only needed by SI3) 849 921 qsr_oce(:,:) = qsr(:,:) 850 922 #endif … … 854 926 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 855 927 CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean 856 CALL iom_put( "qla_oce" , zqla) ! output downward latent heat over the ocean928 CALL iom_put( "qla_oce" , plat ) ! output downward latent heat over the ocean 857 929 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 858 930 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] … … 861 933 ! 862 934 IF ( nn_ice == 0 ) THEN 863 CALL iom_put( "qemp_oce" , qns-zqlw-psen- zqla) ! output downward heat content of E-P over the ocean935 CALL iom_put( "qemp_oce" , qns-zqlw-psen-plat ) ! output downward heat content of E-P over the ocean 864 936 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 865 937 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean … … 869 941 IF(sn_cfctl%l_prtctl) THEN 870 942 CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ') 871 CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') 943 CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen : ' ) 944 CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat : ' ) 945 CALL prt_ctl(tab2d_1=qns , clinfo1=' blk_oce_2: qns : ' ) 872 946 CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ') 873 947 ENDIF … … 883 957 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 884 958 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 885 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag886 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag887 959 !!---------------------------------------------------------------------- 888 960 889 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, p humi, pslp , puice, pvice, ptsui, & ! inputs961 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, pqair, pslp , puice, pvice, ptsui, & ! inputs 890 962 & putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs 891 963 !!--------------------------------------------------------------------- … … 902 974 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s] 903 975 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s] 904 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: p humi! atmospheric wind at T-point [m/s]976 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric wind at T-point [m/s] 905 977 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s] 906 978 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! " … … 915 987 INTEGER :: ji, jj ! dummy loop indices 916 988 REAL(wp) :: zootm_su ! sea-ice surface mean temperature 917 REAL(wp) :: zztmp1, zztmp2 ! temporary arrays 918 REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau) 919 !!--------------------------------------------------------------------- 920 ! 921 989 REAL(wp) :: zztmp1, zztmp2 ! temporary scalars 990 REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! temporary array 991 !!--------------------------------------------------------------------- 992 ! 993 ! LB: ptsui is in K !!! 994 ! 922 995 ! ------------------------------------------------------------ ! 923 996 ! Wind module relative to the moving ice ( U10m - U_ice ) ! … … 925 998 ! C-grid ice dynamics : U & V-points (same as ocean) 926 999 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 927 1000 wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 928 1001 END_2D 929 1002 ! 930 1003 ! Make ice-atm. drag dependent on ice concentration 931 IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations 932 CALL Cdn10_Lupkes2012( Cd_ice ) 933 Ch_ice(:,:) = Cd_ice(:,:) ! momentum and heat transfer coef. are considered identical 934 Ce_ice(:,:) = Cd_ice(:,:) 935 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 936 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 937 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 938 ENDIF 939 940 IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice) 941 IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice) 942 IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice) 943 944 ! local scalars ( place there for vector optimisation purposes) 945 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 1004 1005 1006 SELECT CASE( nblk_ice ) 1007 1008 CASE( np_ice_cst ) 1009 ! Constant bulk transfer coefficients over sea-ice: 1010 Cd_ice(:,:) = rn_Cd_i 1011 Ch_ice(:,:) = rn_Ch_i 1012 Ce_ice(:,:) = rn_Ce_i 1013 ! no height adjustment, keeping zt values: 1014 theta_zu_i(:,:) = ptair(:,:) 1015 q_zu_i(:,:) = pqair(:,:) 1016 1017 CASE( np_ice_an05 ) ! calculate new drag from Lupkes(2015) equations 1018 ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 1019 CALL turb_ice_an05( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, & 1020 & Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 1021 !! 1022 CASE( np_ice_lu12 ) 1023 ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 1024 CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, & 1025 & Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 1026 !! 1027 CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations 1028 ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 1029 CALL turb_ice_lg15( 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 END SELECT 1033 1034 IF( iom_use('Cd_ice').OR.iom_use('Ce_ice').OR.iom_use('Ch_ice').OR.iom_use('taum_ice').OR.iom_use('utau_ice').OR.iom_use('vtau_ice') ) & 1035 & ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice ! 1036 1037 IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp) 1038 IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp) 1039 IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp) 1040 946 1041 947 1042 IF( ln_blk ) THEN … … 950 1045 ! ---------------------------------------------------- ! 951 1046 ! supress moving ice in wind stress computation as we don't know how to do it properly... 952 DO_2D( 0, 1, 0, 1 ) ! at T point 953 putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndi(ji,jj) 954 pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndj(ji,jj) 1047 DO_2D( 0, 1, 0, 1 ) ! at T point 1048 zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj) 1049 putaui(ji,jj) = zztmp1 * pwndi(ji,jj) 1050 pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj) 955 1051 END_2D 1052 1053 !#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!! 1054 IF(iom_use('taum_ice')) CALL iom_put('taum_ice', SQRT( putaui*putaui + pvtaui*pvtaui )*ztmp ) 1055 !#LB: These 2 lines below mostly here for 'STATION_ASF' test-case, otherwize "utau_oi" (U-grid) and vtau_oi" (V-grid) does the job in: [ICE/icedyn_rhg_evp.F90]) 1056 IF(iom_use('utau_ice')) CALL iom_put("utau_ice", putaui*ztmp) ! utau at T-points! 1057 IF(iom_use('vtau_ice')) CALL iom_put("vtau_ice", pvtaui*ztmp) ! vtau at T-points! 1058 956 1059 ! 957 1060 DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean). 958 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 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 959 1063 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 960 1064 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) … … 967 1071 & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) 968 1072 ELSE ! ln_abl 969 zztmp1 = 11637800.0_wp970 zztmp2 = -5897.8_wp971 1073 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 972 pcd_dui(ji,jj) = zcd_dui (ji,jj) 973 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 974 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 975 zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??) 976 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) 977 1077 END_2D 978 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 ) 979 1082 ! 980 1083 IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') … … 983 1086 984 1087 985 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 ) 986 1089 !!--------------------------------------------------------------------- 987 1090 !! *** ROUTINE blk_ice_2 *** … … 999 1102 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 1000 1103 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 1001 REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair 1002 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 1003 1106 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp 1004 REAL(wp), DIMENSION(:,: ), INTENT(in) :: p qlw1107 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pdqlw 1005 1108 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec 1006 1109 REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow 1007 1110 !! 1008 1111 INTEGER :: ji, jj, jl ! dummy loop indices 1009 REAL(wp) :: zst 3! local variable1112 REAL(wp) :: zst, zst3, zsq ! local variable 1010 1113 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 1011 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 1012 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature 1114 REAL(wp) :: zztmp, zzblk, zztmp1, z1_rLsub ! - - 1013 1115 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice 1014 1116 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice … … 1016 1118 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 1017 1119 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3) 1018 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB1019 1120 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2 1020 1121 REAL(wp), DIMENSION(jpi,jpj) :: ztri 1021 1122 !!--------------------------------------------------------------------- 1022 1123 ! 1023 zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars 1024 zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 1025 ! 1026 SELECT CASE( nhumi ) 1027 CASE( np_humi_sph ) 1028 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 1029 CASE( np_humi_dpt ) 1030 zqair(:,:) = q_sat( phumi(:,:), pslp ) 1031 CASE( np_humi_rlh ) 1032 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 1033 END SELECT 1034 ! 1124 zcoef_dqlw = 4._wp * emiss_i * stefan ! local scalars 1125 ! 1126 1035 1127 zztmp = 1. / ( 1. - albo ) 1036 WHERE( ptsu(:,:,:) /= 0._wp ) 1037 z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 1038 ELSEWHERE 1039 z1_st(:,:,:) = 0._wp 1040 END WHERE 1128 dqla_ice(:,:,:) = 0._wp 1129 1041 1130 ! ! ========================== ! 1042 1131 DO jl = 1, jpl ! Loop over ice categories ! 1043 1132 ! ! ========================== ! 1044 DO jj = 1 , jpj 1045 DO ji = 1, jpi 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 1046 1138 ! ----------------------------! 1047 1139 ! I Radiative FLUXES ! 1048 1140 ! ----------------------------! 1049 zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)1050 1141 ! Short Wave (sw) 1051 1142 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 1143 1052 1144 ! Long Wave (lw) 1053 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 1145 zst3 = zst * zst * zst 1146 z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1) 1054 1147 ! lw sensitivity 1055 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst31148 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 1056 1149 1057 1150 ! ----------------------------! … … 1060 1153 1061 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 1062 1159 ! Sensible Heat 1063 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)) 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 1064 1164 ! Latent Heat 1065 zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 1066 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1067 & ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 1068 ! Latent heat sensitivity for ice (Dqla/Dt) 1069 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 1070 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1071 & z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2 1072 ELSE 1073 dqla_ice(ji,jj,jl) = 0._wp 1074 ENDIF 1075 1076 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 1077 z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 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 1078 1173 1079 1174 ! ----------------------------! … … 1083 1178 qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 1084 1179 ! Total non solar heat flux sensitivity for ice 1085 dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) 1086 END DO 1087 ! 1088 END DO 1180 dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ???? 1181 1182 END_2D 1089 1183 ! 1090 1184 END DO … … 1138 1232 ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm 1139 1233 DO jl = 1, jpl 1140 WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1234 WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1141 1235 qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 1142 1236 ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm 1143 1237 qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:) 1144 1238 ELSEWHERE ! zero when hs>0 1145 qtr_ice_top(:,:,jl) = 0._wp 1239 qtr_ice_top(:,:,jl) = 0._wp 1146 1240 END WHERE 1147 1241 ENDDO … … 1182 1276 CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ') 1183 1277 ENDIF 1184 ! 1278 1279 !#LB: 1280 ! air-ice heat flux components that are not written from ice_stp()@icestp.F90: 1281 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 1282 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 1283 IF( iom_use('qlw_ice') ) CALL iom_put( 'qlw_ice', SUM( z_qlw * a_i_b, dim=3 ) ) 1284 !#LB. 1285 1185 1286 END SUBROUTINE blk_ice_2 1186 1287 … … 1278 1379 END SUBROUTINE blk_ice_qcn 1279 1380 1280 1281 SUBROUTINE Cdn10_Lupkes2012( pcd )1282 !!----------------------------------------------------------------------1283 !! *** ROUTINE Cdn10_Lupkes2012 ***1284 !!1285 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m1286 !! to make it dependent on edges at leads, melt ponds and flows.1287 !! After some approximations, this can be resumed to a dependency1288 !! on ice concentration.1289 !!1290 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50)1291 !! with the highest level of approximation: level4, eq.(59)1292 !! The generic drag over a cell partly covered by ice can be re-written as follows:1293 !!1294 !! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu1295 !!1296 !! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59)1297 !! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)1298 !! A is the concentration of ice minus melt ponds (if any)1299 !!1300 !! This new drag has a parabolic shape (as a function of A) starting at1301 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.51302 !! and going down to Cdi(say 1.4e-3) for A=11303 !!1304 !! It is theoretically applicable to all ice conditions (not only MIZ)1305 !! => see Lupkes et al (2013)1306 !!1307 !! ** References : Lupkes et al. JGR 2012 (theory)1308 !! Lupkes et al. GRL 2013 (application to GCM)1309 !!1310 !!----------------------------------------------------------------------1311 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd1312 REAL(wp), PARAMETER :: zCe = 2.23e-03_wp1313 REAL(wp), PARAMETER :: znu = 1._wp1314 REAL(wp), PARAMETER :: zmu = 1._wp1315 REAL(wp), PARAMETER :: zbeta = 1._wp1316 REAL(wp) :: zcoef1317 !!----------------------------------------------------------------------1318 zcoef = znu + 1._wp / ( 10._wp * zbeta )1319 1320 ! generic drag over a cell partly covered by ice1321 !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag1322 !! & Cd_ice * at_i_b(:,:) + & ! pure ice drag1323 !! & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu ! change due to sea-ice morphology1324 1325 ! ice-atm drag1326 pcd(:,:) = rCd_ice + & ! pure ice drag1327 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology1328 1329 END SUBROUTINE Cdn10_Lupkes20121330 1331 1332 SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch )1333 !!----------------------------------------------------------------------1334 !! *** ROUTINE Cdn10_Lupkes2015 ***1335 !!1336 !! ** pUrpose : Alternative turbulent transfert coefficients formulation1337 !! between sea-ice and atmosphere with distinct momentum1338 !! and heat coefficients depending on sea-ice concentration1339 !! and atmospheric stability (no meltponds effect for now).1340 !!1341 !! ** Method : The parameterization is adapted from Lupkes et al. (2015)1342 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,1343 !! it considers specific skin and form drags (Andreas et al. 2010)1344 !! to compute neutral transfert coefficients for both heat and1345 !! momemtum fluxes. Atmospheric stability effect on transfert1346 !! coefficient is also taken into account following Louis (1979).1347 !!1348 !! ** References : Lupkes et al. JGR 2015 (theory)1349 !! Lupkes et al. ECHAM6 documentation 2015 (implementation)1350 !!1351 !!----------------------------------------------------------------------1352 !1353 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K]1354 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa]1355 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient1356 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient1357 REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat1358 !1359 ! ECHAM6 constants1360 REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m]1361 REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m]1362 REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m]1363 REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 411364 REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 411365 REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 131366 REAL(wp), PARAMETER :: zc2 = zc * zc1367 REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 141368 REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 301369 REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 511370 REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 561371 REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 261372 REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 261373 REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma1374 REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp1375 !1376 INTEGER :: ji, jj ! dummy loop indices1377 REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu1378 REAL(wp) :: zrib_o, zrib_i1379 REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice1380 REAL(wp) :: zChn_skin_ice, zChn_form_ice1381 REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw1382 REAL(wp) :: zCdn_form_tmp1383 !!----------------------------------------------------------------------1384 1385 ! Momentum Neutral Transfert Coefficients (should be a constant)1386 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 401387 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 71388 zCdn_ice = zCdn_skin_ice ! Eq. 71389 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32)1390 1391 ! Heat Neutral Transfert Coefficients1392 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. 521393 1394 ! Atmospheric and Surface Variables1395 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin1396 zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg]1397 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg]1398 !1399 DO_2D( 0, 0, 0, 0 )1400 ! Virtual potential temperature [K]1401 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean1402 zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice1403 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu1404 1405 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)1406 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean1407 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice1408 1409 ! Momentum and Heat Neutral Transfert Coefficients1410 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 401411 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 531412 1413 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?)1414 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water1415 z0i = z0_skin_ice ! over ice1416 IF( zrib_o <= 0._wp ) THEN1417 zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 101418 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 261419 & )**zgamma )**z1_gamma1420 ELSE1421 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 121422 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 281423 ENDIF1424 1425 IF( zrib_i <= 0._wp ) THEN1426 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 91427 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 251428 ELSE1429 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 111430 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 271431 ENDIF1432 1433 ! Momentum Transfert Coefficients (Eq. 38)1434 pcd(ji,jj) = zCdn_skin_ice * zfmi + &1435 & 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) )1436 1437 ! Heat Transfert Coefficients (Eq. 49)1438 pch(ji,jj) = zChn_skin_ice * zfhi + &1439 & 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) )1440 !1441 END_2D1442 CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1.0_wp, pch, 'T', 1.0_wp )1443 !1444 END SUBROUTINE Cdn10_Lupkes20151445 1446 1381 #endif 1447 1382
Note: See TracChangeset
for help on using the changeset viewer.
