Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/SBC/sbcblk.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/SBC/sbcblk.F90 (modified) (55 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/SBC/sbcblk.F90
r13501 r14789 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 USE lib_fortran ! to use key_nosignedzero42 USE lib_fortran ! to use key_nosignedzero and glob_sum 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 USE sbcblk_algo_ice_an05 49 USE sbcblk_algo_ice_lu12 50 USE sbcblk_algo_ice_lg15 48 51 #endif 49 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009)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 … … 276 348 ! !- fill the bulk structure with namelist informations 277 349 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 350 sf(jp_wndi )%zsgn = -1._wp ; sf(jp_wndj )%zsgn = -1._wp ! vector field at T point: overwrite default definition of zsgn 351 sf(jp_uoatm)%zsgn = -1._wp ; sf(jp_voatm)%zsgn = -1._wp ! vector field at T point: overwrite default definition of zsgn 352 sf(jp_hpgi )%zsgn = -1._wp ; sf(jp_hpgj )%zsgn = -1._wp ! vector field at T point: overwrite default definition of zsgn 278 353 ! 279 354 DO jfpr= 1, jpfld … … 315 390 END DO 316 391 ! 317 IF( ln_wave ) THEN318 !Activated wave module but neither drag nor stokes drift activated319 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN320 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )321 !drag coefficient read from wave model definable only with mfs bulk formulae and core322 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN323 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')324 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN325 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')326 ENDIF327 ELSE328 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) &329 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', &330 & 'with drag coefficient (ln_cdgw =T) ' , &331 & 'or Stokes Drift (ln_sdw=T) ' , &332 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', &333 & 'or Stokes-Coriolis term (ln_stcori=T)' )334 ENDIF335 !336 392 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient 337 393 rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level … … 341 397 ENDIF 342 398 ! 343 ! set transfer coefficients to default sea-ice values344 Cd_ice(:,:) = rCd_ice345 Ch_ice(:,:) = rCd_ice346 Ce_ice(:,:) = rCd_ice347 399 ! 348 400 IF(lwp) THEN !** Control print … … 350 402 WRITE(numout,*) !* namelist 351 403 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 352 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR404 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 353 405 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 406 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013) ln_COARE_3p6 = ', ln_COARE_3p6 407 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 408 WRITE(numout,*) ' "ANDREAS" algorithm (Andreas et al. 2015) ln_ANDREAS = ', ln_ANDREAS 356 409 WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt 357 410 WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu … … 359 412 WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac 360 413 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 414 WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk 364 415 IF(ln_crt_fbk) THEN … … 374 425 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 375 426 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)' 427 CASE( np_ANDREAS ) ; WRITE(numout,*) ' ==>>> "ANDREAS" algorithm (Andreas et al. 2015)' 376 428 END SELECT 377 429 ! … … 386 438 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 387 439 END SELECT 440 ! 441 !#LB: 442 #if defined key_si3 443 IF( nn_ice > 0 ) THEN 444 WRITE(numout,*) 445 WRITE(numout,*) ' use constant ice-atm bulk transfer coeff. ln_Cx_ice_cst = ', ln_Cx_ice_cst 446 WRITE(numout,*) ' use ice-atm bulk coeff. from Andreas et al., 2005 ln_Cx_ice_AN05 = ', ln_Cx_ice_AN05 447 WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes et al., 2012 ln_Cx_ice_LU12 = ', ln_Cx_ice_LU12 448 WRITE(numout,*) ' use ice-atm bulk coeff. from Lupkes & Gryanik, 2015 ln_Cx_ice_LG15 = ', ln_Cx_ice_LG15 449 ENDIF 450 WRITE(numout,*) 451 SELECT CASE( nblk_ice ) !* Print the choice of bulk algorithm 452 CASE( np_ice_cst ) 453 WRITE(numout,*) ' ==>>> Constant bulk transfer coefficients over sea-ice:' 454 WRITE(numout,*) ' => Cd_ice, Ce_ice, Ch_ice =', REAL(rn_Cd_i,4), REAL(rn_Ce_i,4), REAL(rn_Ch_i,4) 455 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) ) & 456 & CALL ctl_stop( 'Be realistic in your pick of Cd_ice, Ce_ice & Ch_ice ! (0 < Cx < 1.E-2)') 457 CASE( np_ice_an05 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Andreas et al, 2005' 458 CASE( np_ice_lu12 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes et al, 2012' 459 CASE( np_ice_lg15 ) ; WRITE(numout,*) ' ==>>> bulk algo over ice: Lupkes & Gryanik, 2015' 460 END SELECT 461 #endif 462 !#LB. 388 463 ! 389 464 ENDIF … … 428 503 INTEGER, INTENT(in) :: kt ! ocean time step 429 504 !!---------------------------------------------------------------------- 430 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zevp 431 REAL(wp) :: ztmp 505 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zlat, zevp 506 REAL(wp) :: ztst 507 LOGICAL :: llerr 432 508 !!---------------------------------------------------------------------- 433 509 ! … … 436 512 ! Sanity/consistence test on humidity at first time step to detect potential screw-up: 437 513 IF( kt == nit000 ) THEN 438 IF(lwp) WRITE(numout,*) '' 439 #if defined key_agrif 440 IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 441 #else 442 ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain 443 IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points! 444 ztmp = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztmp ! mean humidity over ocean on proc 445 SELECT CASE( nhumi ) 446 CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!) 447 IF( (ztmp < 0._wp) .OR. (ztmp > 0.065) ) ztmp = -1._wp 448 CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K] 449 IF( (ztmp < 110._wp).OR.(ztmp > 320._wp) ) ztmp = -1._wp 450 CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%] 451 IF( (ztmp < 0._wp) .OR.(ztmp > 100._wp) ) ztmp = -1._wp 452 END SELECT 453 IF(ztmp < 0._wp) THEN 454 IF (lwp) WRITE(numout,'(" Mean humidity value found on proc #",i6.6," is: ",f10.5)') narea, ztmp 455 CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', & 456 & ' ==> check the unit in your input files' , & 457 & ' ==> check consistence of namelist choice: specific? relative? dew-point?', & 458 & ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' ) 459 END IF 460 END IF 461 IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 462 #endif 463 IF(lwp) WRITE(numout,*) '' 464 END IF !IF( kt == nit000 ) 514 ! mean humidity over ocean on proc 515 ztst = glob_sum( 'sbcblk', sf(jp_humi)%fnow(:,:,1) * e1e2t(:,:) * tmask(:,:,1) ) / glob_sum( 'sbcblk', e1e2t(:,:) * tmask(:,:,1) ) 516 llerr = .FALSE. 517 SELECT CASE( nhumi ) 518 CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!) 519 IF( (ztst < 0._wp) .OR. (ztst > 0.065_wp) ) llerr = .TRUE. 520 CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K] 521 IF( (ztst < 110._wp) .OR. (ztst > 320._wp) ) llerr = .TRUE. 522 CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%] 523 IF( (ztst < 0._wp) .OR. (ztst > 100._wp) ) llerr = .TRUE. 524 END SELECT 525 IF(llerr) THEN 526 WRITE(ctmp1,'(" Error on mean humidity value: ",f10.5)') ztst 527 CALL ctl_stop( 'STOP', ctmp1, 'Something is wrong with air humidity!!!', & 528 & ' ==> check the unit in your input files' , & 529 & ' ==> check consistence of namelist choice: specific? relative? dew-point?', & 530 & ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' ) 531 ENDIF 532 IF(lwp) THEN 533 WRITE(numout,*) '' 534 WRITE(numout,*) ' Global mean humidity at kt = nit000: ', ztst 535 WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 536 WRITE(numout,*) '' 537 ENDIF 538 ENDIF !IF( kt == nit000 ) 465 539 ! ! compute the surface ocean fluxes using bulk formulea 466 540 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 541 542 ! Specific humidity of air at z=rn_zqt ! 543 SELECT CASE( nhumi ) 544 CASE( np_humi_sph ) 545 q_air_zt(:,:) = sf(jp_humi )%fnow(:,:,1) ! what we read in file is already a spec. humidity! 546 CASE( np_humi_dpt ) 547 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing q_air out of dew-point and P !' 548 q_air_zt(:,:) = q_sat( sf(jp_humi )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) 549 CASE( np_humi_rlh ) 550 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => computing q_air out of RH, t_air and slp !' !LBrm 551 q_air_zt(:,:) = q_air_rh( 0.01_wp*sf(jp_humi )%fnow(:,:,1), & 552 & sf(jp_tair )%fnow(:,:,1), sf(jp_slp )%fnow(:,:,1) ) !#LB: 0.01 => RH is % percent in file 553 END SELECT 554 555 ! POTENTIAL temperature of air at z=rn_zqt 556 ! based on adiabatic lapse-rate (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 557 ! (most reanalysis products provide absolute temp., not potential temp.) 558 IF( ln_tair_pot ) THEN 559 ! temperature read into file is already potential temperature, do nothing... 560 theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) 561 ELSE 562 ! temperature read into file is ABSOLUTE temperature (that's the case for ECMWF products for example...) 563 IF((kt==nit000).AND.lwp) WRITE(numout,*) ' *** sbc_blk() => air temperature converted from ABSOLUTE to POTENTIAL!' 564 theta_air_zt(:,:) = sf(jp_tair )%fnow(:,:,1) + gamma_moist( sf(jp_tair )%fnow(:,:,1), q_air_zt(:,:) ) * rn_zqt 565 ENDIF 566 ! 467 567 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),& ! <<= in568 & theta_air_zt(:,:), q_air_zt(:,:), & ! <<= in 469 569 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 470 570 & sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in 471 571 & 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),& ! <<= in572 & tsk_m, zssq, zcd_du, zsen, zlat, zevp ) ! =>> out 573 574 CALL blk_oce_2( theta_air_zt(:,:), & ! <<= in 475 575 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1), & ! <<= in 476 576 & sf(jp_snow )%fnow(:,:,1), tsk_m, & ! <<= in 477 & zsen, z evp )! <=> in out577 & zsen, zlat, zevp ) ! <=> in out 478 578 ENDIF 479 579 ! … … 486 586 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 487 587 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 588 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) !#LB: should it be POTENTIAL temperature instead ???? 589 !tatm_ice(:,:) = theta_air_zt(:,:) !#LB: THIS! ? 590 591 qatm_ice(:,:) = q_air_zt(:,:) !#LB: 498 592 499 593 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac … … 507 601 508 602 509 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, p humi, & ! inp603 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, pqair, & ! inp 510 604 & pslp , pst , pu , pv, & ! inp 511 & puatm, pvatm, p qsr , pqlw ,& ! inp512 & ptsk , pssq , pcd_du, psen, p evp )! out605 & puatm, pvatm, pdqsr , pdqlw , & ! inp 606 & ptsk , pssq , pcd_du, psen, plat, pevp ) ! out 513 607 !!--------------------------------------------------------------------- 514 608 !! *** ROUTINE blk_oce_1 *** … … 523 617 !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) 524 618 !! - 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) 619 !! - psen : sensible heat flux (W/m^2) 620 !! - plat : latent heat flux (W/m^2) 621 !! - pevp : evaporation (mm/s) #lolo 527 622 !!--------------------------------------------------------------------- 528 623 INTEGER , INTENT(in ) :: kt ! time step index 529 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s]530 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]624 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at T-point [m/s] 625 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at T-point [m/s] 626 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqair ! specific humidity at T-points [kg/kg] 532 627 REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin] 533 628 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa] … … 537 632 REAL(wp), INTENT(in ), DIMENSION(:,:) :: puatm ! surface current seen by the atm at T-point (i-component) [m/s] 538 633 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 !634 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqsr ! downwelling solar (shortwave) radiation at surface [W/m^2] 635 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2] 541 636 REAL(wp), INTENT( out), DIMENSION(:,:) :: ptsk ! skin temp. (or SST if CS & WL not used) [Celsius] 542 637 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] 638 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du 639 REAL(wp), INTENT( out), DIMENSION(:,:) :: psen 640 REAL(wp), INTENT( out), DIMENSION(:,:) :: plat 641 REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp 546 642 ! 547 643 INTEGER :: ji, jj ! dummy loop indices … … 553 649 REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point 554 650 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 651 REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean 558 652 REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean 559 653 REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean 560 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux561 654 REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2 562 655 !!--------------------------------------------------------------------- … … 597 690 zztmp = 1. - albo 598 691 IF( ln_dm2dc ) THEN 599 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)692 qsr(:,:) = zztmp * sbc_dcy( pdqsr(:,:) ) * tmask(:,:,1) 600 693 ELSE 601 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)694 qsr(:,:) = zztmp * pdqsr(:,:) * tmask(:,:,1) 602 695 ENDIF 603 696 … … 616 709 ENDIF 617 710 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 711 !! Time to call the user-selected bulk parameterization for 649 712 !! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more... … … 651 714 652 715 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 716 CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 717 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , & 718 & nb_iter=nn_iter_algo ) 719 ! 656 720 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 721 CALL turb_coare3p0( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 722 & ln_skin_cs, ln_skin_wl, & 723 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 724 & nb_iter=nn_iter_algo, & 725 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 726 ! 661 727 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 728 CALL turb_coare3p6( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 729 & ln_skin_cs, ln_skin_wl, & 730 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 731 & nb_iter=nn_iter_algo, & 732 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 733 ! 666 734 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 735 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 736 & ln_skin_cs, ln_skin_wl, & 737 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu, & 738 & nb_iter=nn_iter_algo, & 739 & Qsw=qsr(:,:), rad_lw=pdqlw(:,:), slp=pslp(:,:) ) 740 ! 741 CASE( np_ANDREAS ) 742 CALL turb_andreas ( rn_zqt, rn_zu, ptsk, ptair, pssq, pqair, wndm, & 743 & zcd_oce, zch_oce, zce_oce, theta_zu, q_zu, zU_zu , & 744 & nb_iter=nn_iter_algo ) 745 ! 671 746 CASE DEFAULT 672 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formulaselected' )673 747 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk parameterizaton selected' ) 748 ! 674 749 END SELECT 675 750 676 751 IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1)) 677 752 IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1)) 678 753 IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1)) 679 754 !! 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=zu755 IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ptair-rt0) * tmask(:,:,1)) ! potential temperature at z=zt 756 IF( iom_use('q_zt') ) CALL iom_put("q_zt", pqair * tmask(:,:,1)) ! specific humidity " 757 IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (theta_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu 683 758 IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity " 684 759 IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0 685 760 IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu 686 761 687 762 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 688 763 !! ptsk and pssq have been updated!!! … … 696 771 END IF 697 772 698 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbc blk_phy.F90773 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbc_phy.F90 699 774 ! ------------------------------------------------------------- 700 775 … … 706 781 psen(ji,jj) = zztmp * zch_oce(ji,jj) 707 782 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 708 rhoa(ji,jj) = rho_air( ptair(ji,jj), p humi(ji,jj), pslp(ji,jj) )783 rhoa(ji,jj) = rho_air( ptair(ji,jj), pqair(ji,jj), pslp(ji,jj) ) 709 784 END_2D 710 785 ELSE !== BLK formulation ==! turbulent fluxes computation 711 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t _zu(:,:), q_zu(:,:), &786 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), & 712 787 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 713 788 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 714 & taum(:,:), psen(:,:), zqla(:,:), &789 & taum(:,:), psen(:,:), plat(:,:), & 715 790 & pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac ) 716 791 717 zqla(:,:) = zqla(:,:) * tmask(:,:,1)718 792 psen(:,:) = psen(:,:) * tmask(:,:,1) 793 plat(:,:) = plat(:,:) * tmask(:,:,1) 719 794 taum(:,:) = taum(:,:) * tmask(:,:,1) 720 795 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 721 796 722 797 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 723 724 798 IF( wndm(ji,jj) > 0._wp ) THEN 799 zztmp = taum(ji,jj) / wndm(ji,jj) 725 800 #if defined key_cyclone 726 727 801 ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj) 802 ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj) 728 803 #else 729 730 804 ztau_i(ji,jj) = zztmp * pwndi(ji,jj) 805 ztau_j(ji,jj) = zztmp * pwndj(ji,jj) 731 806 #endif 732 733 734 ztau_j(ji,jj) = 0._wp735 807 ELSE 808 ztau_i(ji,jj) = 0._wp 809 ztau_j(ji,jj) = 0._wp 810 ENDIF 736 811 END_2D 737 812 … … 757 832 758 833 IF( ln_crt_fbk ) THEN 759 CALL lbc_lnk _multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1.)834 CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp, taum, 'T', 1._wp ) 760 835 ELSE 761 CALL lbc_lnk _multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1.)836 CALL lbc_lnk( 'sbcblk', utau, 'U', -1._wp, vtau, 'V', -1._wp ) 762 837 ENDIF 763 838 764 CALL iom_put( "taum_oce", taum ) ! output wind stress module 839 ! Saving open-ocean wind-stress (module and components) on T-points: 840 CALL iom_put( "taum_oce", taum(:,:)*tmask(:,:,1) ) ! output wind stress module 841 !#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]) 842 CALL iom_put( "utau_oce", ztau_i(:,:)*tmask(:,:,1) ) ! utau at T-points! 843 CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) ) ! vtau at T-points! 765 844 766 845 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 ) 846 CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ') 847 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') 848 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & 849 & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 850 CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd : ') 770 851 ENDIF 771 852 ! 772 853 ENDIF !IF( ln_abl ) 773 854 774 855 ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius 775 856 776 857 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 777 858 CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius 778 859 CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference... 779 860 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 861 ! 787 862 END SUBROUTINE blk_oce_1 788 863 789 864 790 SUBROUTINE blk_oce_2( ptair, p qsr, pqlw, pprec,& ! <<= in791 & psnow, ptsk, psen, pevp ) ! <<= in865 SUBROUTINE blk_oce_2( ptair, pdqlw, pprec, psnow, & ! <<= in 866 & ptsk, psen, plat, pevp ) ! <<= in 792 867 !!--------------------------------------------------------------------- 793 868 !! *** ROUTINE blk_oce_2 *** … … 805 880 !! - emp : evaporation minus precipitation (kg/m2/s) 806 881 !!--------------------------------------------------------------------- 807 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair 808 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr 809 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw 882 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair ! potential temperature of air #LB: confirm! 883 REAL(wp), INTENT(in), DIMENSION(:,:) :: pdqlw ! downwelling longwave radiation at surface [W/m^2] 810 884 REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec 811 885 REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow 812 886 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptsk ! SKIN surface temperature [Celsius] 813 887 REAL(wp), INTENT(in), DIMENSION(:,:) :: psen 888 REAL(wp), INTENT(in), DIMENSION(:,:) :: plat 814 889 REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp 815 890 ! 816 891 INTEGER :: ji, jj ! dummy loop indices 817 892 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 821 !!--------------------------------------------------------------------- 822 ! 823 ! local scalars ( place there for vector optimisation purposes) 824 825 826 ztskk(:,:) = ptsk(:,:) + rt0 ! => ptsk in Kelvin rather than Celsius 827 893 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! net long wave radiative heat flux 894 REAL(wp), DIMENSION(jpi,jpj) :: zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg) 895 !!--------------------------------------------------------------------- 896 ! 897 ! Heat content per unit mass (J/kg) 898 zcptrain(:,:) = ( ptair - rt0 ) * rcp * tmask(:,:,1) 899 zcptsnw (:,:) = ( MIN( ptair, rt0 ) - rt0 ) * rcpi * tmask(:,:,1) 900 zcptn (:,:) = ptsk * rcp * tmask(:,:,1) 901 ! 828 902 ! ----------------------------------------------------------------------------- ! 829 903 ! III Net longwave radiative FLUX ! 830 904 ! ----------------------------------------------------------------------------- ! 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 905 !! #LB: now moved after Turbulent fluxes because must use the skin temperature rather than bulk SST 906 !! (ptsk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 907 zqlw(:,:) = qlw_net( pdqlw(:,:), ptsk(:,:)+rt0 ) 849 908 850 909 ! ----------------------------------------------------------------------------- ! … … 852 911 ! ----------------------------------------------------------------------------- ! 853 912 ! 854 emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.) 855 & - pprec(:,:) * rn_pfac ) * tmask(:,:,1) 856 ! 857 qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar 858 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 859 & - pevp(:,:) * ptsk(:,:) * rcp & ! remove evap heat content at SST 860 & + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair 861 & * ( ptair(:,:) - rt0 ) * rcp & 862 & + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 863 & * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 913 emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac ) * tmask(:,:,1) ! mass flux (evap. - precip.) 914 ! 915 qns(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) & ! Downward Non Solar 916 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 917 & - pevp(:,:) * zcptn(:,:) & ! remove evap heat content at SST 918 & + ( pprec(:,:) - psnow(:,:) ) * rn_pfac * zcptrain(:,:) & ! add liquid precip heat content at Tair 919 & + psnow(:,:) * rn_pfac * zcptsnw(:,:) ! add solid precip heat content at min(Tair,Tsnow) 864 920 qns(:,:) = qns(:,:) * tmask(:,:,1) 865 921 ! 866 922 #if defined key_si3 867 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3)923 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + plat(:,:) ! non solar without emp (only needed by SI3) 868 924 qsr_oce(:,:) = qsr(:,:) 869 925 #endif … … 873 929 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 874 930 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 ocean931 CALL iom_put( "qla_oce" , plat ) ! output downward latent heat over the ocean 876 932 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 877 933 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] … … 880 936 ! 881 937 IF ( nn_ice == 0 ) THEN 882 CALL iom_put( "qemp_oce" , qns-zqlw-psen- zqla) ! output downward heat content of E-P over the ocean938 CALL iom_put( "qemp_oce" , qns-zqlw-psen-plat ) ! output downward heat content of E-P over the ocean 883 939 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 884 940 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean … … 888 944 IF(sn_cfctl%l_prtctl) THEN 889 945 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 : ') 946 CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen : ' ) 947 CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat : ' ) 948 CALL prt_ctl(tab2d_1=qns , clinfo1=' blk_oce_2: qns : ' ) 891 949 CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ') 892 950 ENDIF … … 902 960 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 903 961 !! 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 962 !!---------------------------------------------------------------------- 907 963 908 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, p humi, pslp , puice, pvice, ptsui, & ! inputs964 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, pqair, pslp , puice, pvice, ptsui, & ! inputs 909 965 & putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs 910 966 !!--------------------------------------------------------------------- … … 921 977 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s] 922 978 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]979 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pqair ! atmospheric wind at T-point [m/s] 924 980 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s] 925 981 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! " … … 934 990 INTEGER :: ji, jj ! dummy loop indices 935 991 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 992 REAL(wp) :: zztmp1, zztmp2 ! temporary scalars 993 REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! temporary array 994 !!--------------------------------------------------------------------- 995 ! 996 ! LB: ptsui is in K !!! 997 ! 941 998 ! ------------------------------------------------------------ ! 942 999 ! Wind module relative to the moving ice ( U10m - U_ice ) ! … … 948 1005 ! 949 1006 ! 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(:,:) 1007 1008 1009 SELECT CASE( nblk_ice ) 1010 1011 CASE( np_ice_cst ) 1012 ! Constant bulk transfer coefficients over sea-ice: 1013 Cd_ice(:,:) = rn_Cd_i 1014 Ch_ice(:,:) = rn_Ch_i 1015 Ce_ice(:,:) = rn_Ce_i 1016 ! no height adjustment, keeping zt values: 1017 theta_zu_i(:,:) = ptair(:,:) 1018 q_zu_i(:,:) = pqair(:,:) 1019 1020 CASE( np_ice_an05 ) ! calculate new drag from Lupkes(2015) equations 1021 ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 1022 CALL turb_ice_an05( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, & 1023 & Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 1024 !! 1025 CASE( np_ice_lu12 ) 1026 ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 1027 CALL turb_ice_lu12( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, & 1028 & Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 1029 !! 1030 CASE( np_ice_lg15 ) ! calculate new drag from Lupkes(2015) equations 1031 ztmp(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ! temporary array for SSQ 1032 CALL turb_ice_lg15( rn_zqt, rn_zu, ptsui, ptair, ztmp, pqair, wndm_ice, fr_i, & 1033 & Cd_ice, Ch_ice, Ce_ice, theta_zu_i, q_zu_i ) 1034 !! 1035 END SELECT 1036 1037 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') ) & 1038 & ztmp(:,:) = ( 1._wp - MAX(0._wp, SIGN( 1._wp, 1.E-6_wp - fr_i )) )*tmask(:,:,1) ! mask for presence of ice ! 1039 1040 IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice*ztmp) 1041 IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice*ztmp) 1042 IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice*ztmp) 1043 965 1044 966 1045 IF( ln_blk ) THEN … … 969 1048 ! ---------------------------------------------------- ! 970 1049 ! 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) 1050 DO_2D( 0, 1, 0, 1 ) ! at T point 1051 zztmp1 = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj) 1052 putaui(ji,jj) = zztmp1 * pwndi(ji,jj) 1053 pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj) 974 1054 END_2D 1055 1056 !#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!! 1057 IF(iom_use('taum_ice')) CALL iom_put('taum_ice', SQRT( putaui*putaui + pvtaui*pvtaui )*ztmp ) 1058 !#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]) 1059 IF(iom_use('utau_ice')) CALL iom_put("utau_ice", putaui*ztmp) ! utau at T-points! 1060 IF(iom_use('vtau_ice')) CALL iom_put("vtau_ice", pvtaui*ztmp) ! vtau at T-points! 1061 975 1062 ! 976 1063 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 1064 !#LB: QUESTION?? so SI3 expects wind stress vector to be provided at U & V points? Not at T-points ? 1065 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 978 1066 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 979 1067 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) … … 981 1069 pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) ) 982 1070 END_2D 983 CALL lbc_lnk _multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )1071 CALL lbc_lnk( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp ) 984 1072 ! 985 1073 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' & 986 1074 & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) 987 1075 ELSE ! ln_abl 988 zztmp1 = 11637800.0_wp989 zztmp2 = -5897.8_wp990 1076 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) 1077 pcd_dui(ji,jj) = wndm_ice(ji,jj) * Cd_ice(ji,jj) 1078 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 1079 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 996 1080 END_2D 997 ENDIF 1081 !#LB: 1082 pssqi(:,:) = q_sat( ptsui(:,:), pslp(:,:), l_ice=.TRUE. ) ; ! more accurate way to obtain ssq ! 1083 !#LB. 1084 ENDIF !IF( ln_blk ) 998 1085 ! 999 1086 IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') … … 1002 1089 1003 1090 1004 SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, p humi, pslp, pqlw, pprec, psnow )1091 SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, pqair, pslp, pdqlw, pprec, psnow ) 1005 1092 !!--------------------------------------------------------------------- 1006 1093 !! *** ROUTINE blk_ice_2 *** … … 1018 1105 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 1019 1106 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 1020 REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair 1021 REAL(wp), DIMENSION(:,: ), INTENT(in) :: p humi1107 REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair ! potential temperature of air #LB: okay ??? 1108 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqair ! specific humidity of air 1022 1109 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp 1023 REAL(wp), DIMENSION(:,: ), INTENT(in) :: p qlw1110 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pdqlw 1024 1111 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec 1025 1112 REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow 1026 1113 !! 1027 1114 INTEGER :: ji, jj, jl ! dummy loop indices 1028 REAL(wp) :: zst 3! local variable1115 REAL(wp) :: zst, zst3, zsq ! local variable 1029 1116 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 1117 REAL(wp) :: zztmp, zzblk, zztmp1, z1_rLsub ! - - 1032 1118 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice 1033 1119 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice … … 1035 1121 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 1036 1122 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 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp21039 1123 REAL(wp), DIMENSION(jpi,jpj) :: ztri 1040 !!--------------------------------------------------------------------- 1041 ! 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 REAL(wp), DIMENSION(jpi,jpj) :: zcptrain, zcptsnw, zcptn ! Heat content per unit mass (J/kg) 1125 !!--------------------------------------------------------------------- 1126 ! 1127 zcoef_dqlw = 4._wp * emiss_i * stefan ! local scalars 1128 ! 1129 1054 1130 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 1131 dqla_ice(:,:,:) = 0._wp 1132 1133 ! Heat content per unit mass (J/kg) 1134 zcptrain(:,:) = ( ptair - rt0 ) * rcp * tmask(:,:,1) 1135 zcptsnw (:,:) = ( MIN( ptair, rt0 ) - rt0 ) * rcpi * tmask(:,:,1) 1136 zcptn (:,:) = sst_m * rcp * tmask(:,:,1) 1137 ! 1060 1138 ! ! ========================== ! 1061 1139 DO jl = 1, jpl ! Loop over ice categories ! 1062 1140 ! ! ========================== ! 1063 DO jj = 1 , jpj 1064 DO ji = 1, jpi 1141 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1142 1143 zst = ptsu(ji,jj,jl) ! surface temperature of sea-ice [K] 1144 zsq = q_sat( zst, pslp(ji,jj), l_ice=.TRUE. ) ! surface saturation specific humidity when ice present 1145 1065 1146 ! ----------------------------! 1066 1147 ! I Radiative FLUXES ! 1067 1148 ! ----------------------------! 1068 zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)1069 1149 ! Short Wave (sw) 1070 1150 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 1151 1071 1152 ! Long Wave (lw) 1072 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 1153 zst3 = zst * zst * zst 1154 z_qlw(ji,jj,jl) = emiss_i * ( pdqlw(ji,jj) - stefan * zst * zst3 ) * tmask(ji,jj,1) 1073 1155 ! lw sensitivity 1074 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst31156 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 1075 1157 1076 1158 ! ----------------------------! … … 1079 1161 1080 1162 ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 1163 1164 ! Common term in bulk F. equations... 1165 zzblk = rhoa(ji,jj) * wndm_ice(ji,jj) 1166 1081 1167 ! 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)) 1168 zztmp1 = zzblk * rCp_air * Ch_ice(ji,jj) 1169 z_qsb (ji,jj,jl) = zztmp1 * (zst - theta_zu_i(ji,jj)) 1170 z_dqsb(ji,jj,jl) = zztmp1 ! ==> Qsens sensitivity (Dqsb_ice/Dtn_ice) 1171 1083 1172 ! 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) 1173 zztmp1 = zzblk * rLsub * Ce_ice(ji,jj) 1174 qla_ice(ji,jj,jl) = MAX( zztmp1 * (zsq - q_zu_i(ji,jj)) , 0._wp ) ! #LB: only sublimation (and not condensation) ??? 1175 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) 1176 ! !#LB: dq_sat_dt_ice() in "sbc_phy.F90" 1177 !#LB: without this unjustified "condensation sensure": 1178 !qla_ice( ji,jj,jl) = zztmp1 * (zsq - q_zu_i(ji,jj)) 1179 !dqla_ice(ji,jj,jl) = zztmp1 * dq_sat_dt_ice(zst, pslp(ji,jj)) ! ==> Qlat sensitivity (dQlat/dT) 1180 1097 1181 1098 1182 ! ----------------------------! … … 1102 1186 qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 1103 1187 ! 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 1188 dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) !#LB: correct signs ???? 1189 1190 END_2D 1108 1191 ! 1109 1192 END DO … … 1128 1211 1129 1212 ! --- heat flux associated with emp --- ! 1130 qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 1131 & + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair 1132 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) 1133 & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 1134 qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only) 1135 & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 1213 qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * zcptn(:,:) & ! evap at sst 1214 & + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) & ! liquid precip at Tair 1215 & + sprecip(:,:) * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - rLfus ) ! solid precip at min(Tair,Tsnow) 1216 qemp_ice(:,:) = sprecip(:,:) * zsnw * ( zcptsnw (:,:) - rLfus ) ! solid precip (only) 1136 1217 1137 1218 ! --- total solar and non solar fluxes --- ! … … 1141 1222 1142 1223 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 1143 qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1224 qprec_ice(:,:) = rhos * ( zcptsnw(:,:) - rLfus ) 1144 1225 1145 1226 ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- … … 1157 1238 ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm 1158 1239 DO jl = 1, jpl 1159 WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1240 WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1160 1241 qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 1161 1242 ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm 1162 1243 qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:) 1163 1244 ELSEWHERE ! zero when hs>0 1164 qtr_ice_top(:,:,jl) = 0._wp 1245 qtr_ice_top(:,:,jl) = 0._wp 1165 1246 END WHERE 1166 1247 ENDDO … … 1173 1254 ! 1174 1255 IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 1175 ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 1176 IF( iom_use('evap_ao_cea' ) ) CALL iom_put( 'evap_ao_cea' , ztmp(:,:) * tmask(:,:,1) ) ! ice-free oce evap (cell average) 1177 IF( iom_use('hflx_evap_cea') ) CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) ) ! heat flux from evap (cell average) 1178 ENDIF 1179 IF( iom_use('hflx_rain_cea') ) THEN 1180 ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) ) 1181 IF( iom_use('hflx_rain_cea') ) CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) ) ! heat flux from rain (cell average) 1182 ENDIF 1183 IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea') ) THEN 1184 WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) 1185 ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 1186 ELSEWHERE 1187 ztmp(:,:) = rcp * sst_m(:,:) 1188 ENDWHERE 1189 ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 1190 IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea' , ztmp2(:,:) ) ! heat flux from snow (cell average) 1191 IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 1192 IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) * zsnw(:,:) ) ! heat flux from snow (over ice) 1256 CALL iom_put( 'evap_ao_cea' , zevap(:,:) * ( 1._wp - at_i_b(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average) 1257 CALL iom_put( 'hflx_evap_cea', zevap(:,:) * ( 1._wp - at_i_b(:,:) ) * tmask(:,:,1) * zcptn(:,:) ) ! heat flux from evap (cell average) 1258 ENDIF 1259 IF( iom_use('rain') .OR. iom_use('rain_ao_cea') .OR. iom_use('hflx_rain_cea') ) THEN 1260 CALL iom_put( 'rain' , tprecip(:,:) - sprecip(:,:) ) ! liquid precipitation 1261 CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * ( 1._wp - at_i_b(:,:) ) ) ! liquid precipitation over ocean (cell average) 1262 CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) ) ! heat flux from rain (cell average) 1263 ENDIF 1264 IF( iom_use('snow_ao_cea') .OR. iom_use('snow_ai_cea') .OR. & 1265 & iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea') ) THEN 1266 CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) ) ) ! Snow over ice-free ocean (cell average) 1267 CALL iom_put( 'snow_ai_cea' , sprecip(:,:) * zsnw(:,:) ) ! Snow over sea-ice (cell average) 1268 CALL iom_put( 'hflx_snow_cea' , sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) ) ! heat flux from snow (cell average) 1269 CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 1270 CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) * zsnw(:,:) ) ! heat flux from snow (over ice) 1271 ENDIF 1272 IF( iom_use('hflx_prec_cea') ) THEN ! heat flux from precip (cell average) 1273 CALL iom_put('hflx_prec_cea' , sprecip(:,:) * ( zcptsnw (:,:) - rLfus ) & 1274 & + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) ) 1275 ENDIF 1276 IF( iom_use('subl_ai_cea') .OR. iom_use('hflx_subl_cea') ) THEN 1277 CALL iom_put( 'subl_ai_cea' , SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) * tmask(:,:,1) ) ! Sublimation over sea-ice (cell average) 1278 CALL iom_put( 'hflx_subl_cea', SUM( a_i_b(:,:,:) * qevap_ice(:,:,:), dim=3 ) * tmask(:,:,1) ) ! Heat flux from sublimation (cell average) 1193 1279 ENDIF 1194 1280 ! … … 1201 1287 CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ') 1202 1288 ENDIF 1203 ! 1289 1290 !#LB: 1291 ! air-ice heat flux components that are not written from ice_stp()@icestp.F90: 1292 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 1293 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 1294 IF( iom_use('qlw_ice') ) CALL iom_put( 'qlw_ice', SUM( z_qlw * a_i_b, dim=3 ) ) 1295 !#LB. 1296 1204 1297 END SUBROUTINE blk_ice_2 1205 1298 … … 1297 1390 END SUBROUTINE blk_ice_qcn 1298 1391 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 1465 1392 #endif 1466 1393
Note: See TracChangeset
for help on using the changeset viewer.
