- Timestamp:
- 2020-09-29T12:41:06+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/r12377_ticket2386
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/SBC/sbcblk.F90 (modified) (30 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/r12377_ticket2386
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 8 9 9 # SETTE 10 ^/utils/CI/sette@ HEADsette10 ^/utils/CI/sette@13507 sette
-
- Property svn:externals
-
NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk.F90
r12511 r13540 44 44 USE lib_fortran ! to use key_nosignedzero 45 45 #if defined key_si3 46 USE ice , ONLY : jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif47 USE ice thd_dh ! for CALL ice_thd_snwblow46 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 USE icevar ! for CALL ice_var_snwblow 48 48 #endif 49 49 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) … … 74 74 #endif 75 75 76 INTEGER , PUBLIC :: jpfld ! maximum number of files to read 77 INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point 78 INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point 79 INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin) 80 INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 81 INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2) 82 INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2) 83 INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 84 INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 85 INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa) 86 INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 87 INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 88 76 INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point 77 INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point 78 INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin) 79 INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 80 INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2) 81 INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2) 82 INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 83 INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 84 INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa) 85 INTEGER , PUBLIC, PARAMETER :: jp_uoatm = 10 ! index of surface current (i-component) 86 ! ! seen by the atmospheric forcing (m/s) at T-point 87 INTEGER , PUBLIC, PARAMETER :: jp_voatm = 11 ! index of surface current (j-component) 88 ! ! seen by the atmospheric forcing (m/s) at T-point 89 INTEGER , PUBLIC, PARAMETER :: jp_cc = 12 ! index of cloud cover (-) range:0-1 90 INTEGER , PUBLIC, PARAMETER :: jp_hpgi = 13 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 91 INTEGER , PUBLIC, PARAMETER :: jp_hpgj = 14 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 92 INTEGER , PUBLIC, PARAMETER :: jpfld = 14 ! maximum number of files to read 93 94 ! Warning: keep this structure allocatable for Agrif... 89 95 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read) 90 96 … … 98 104 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 99 105 ! 106 LOGICAL :: ln_crt_fbk ! Add surface current feedback to the wind stress computation (Renault et al. 2020) 107 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 ! 109 ! 100 110 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 101 111 REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation 102 REAL(wp), PUBLIC :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress103 112 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 104 113 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 105 114 ! 106 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice107 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme)108 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme)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) 109 118 110 119 LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB … … 113 122 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 114 123 LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB 124 LOGICAL :: ln_tpot !!GS: flag to compute or not potential temperature 115 125 ! 116 126 INTEGER :: nhumi ! choice of the bulk algorithm … … 162 172 !! 163 173 CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files 164 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 165 TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read 166 TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " " 167 TYPE(FLD_N) :: sn_slp , sn_hpgi, sn_hpgj ! " " 174 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 175 TYPE(FLD_N) :: sn_wndi, sn_wndj , sn_humi, sn_qsr ! informations about the fields to be read 176 TYPE(FLD_N) :: sn_qlw , sn_tair , sn_prec, sn_snow ! " " 177 TYPE(FLD_N) :: sn_slp , sn_uoatm, sn_voatm ! " " 178 TYPE(FLD_N) :: sn_cc, sn_hpgi, sn_hpgj ! " " 179 INTEGER :: ipka ! number of levels in the atmospheric variable 168 180 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 169 & sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj, & 181 & sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm, & 182 & sn_cc, sn_hpgi, sn_hpgj, & 170 183 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 171 184 & cn_dir , rn_zqt, rn_zu, & 172 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, & 185 & rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15, ln_tpot, & 186 & ln_crt_fbk, rn_stau_a, rn_stau_b, & ! current feedback 173 187 & ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB 174 188 !!--------------------------------------------------------------------- … … 242 256 ! !* set the bulk structure 243 257 ! !- store namelist information in an array 244 IF( ln_blk ) jpfld = 9 245 IF( ln_abl ) jpfld = 11 246 ALLOCATE( slf_i(jpfld) ) 247 ! 248 slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj 249 slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw 250 slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi 251 slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow 252 slf_i(jp_slp ) = sn_slp 253 IF( ln_abl ) THEN 254 slf_i(jp_hpgi) = sn_hpgi ; slf_i(jp_hpgj) = sn_hpgj 255 END IF 258 ! 259 slf_i(jp_wndi ) = sn_wndi ; slf_i(jp_wndj ) = sn_wndj 260 slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw 261 slf_i(jp_tair ) = sn_tair ; slf_i(jp_humi ) = sn_humi 262 slf_i(jp_prec ) = sn_prec ; slf_i(jp_snow ) = sn_snow 263 slf_i(jp_slp ) = sn_slp ; slf_i(jp_cc ) = sn_cc 264 slf_i(jp_uoatm) = sn_uoatm ; slf_i(jp_voatm) = sn_voatm 265 slf_i(jp_hpgi ) = sn_hpgi ; slf_i(jp_hpgj ) = sn_hpgj 266 ! 267 IF( .NOT. ln_abl ) THEN ! force to not use jp_hpgi and jp_hpgj, should already be done in namelist_* but we never know... 268 slf_i(jp_hpgi)%clname = 'NOT USED' 269 slf_i(jp_hpgj)%clname = 'NOT USED' 270 ENDIF 256 271 ! 257 272 ! !- allocate the bulk structure … … 264 279 DO jfpr= 1, jpfld 265 280 ! 266 IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to zero) 267 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 268 sf(jfpr)%fnow(:,:,1) = 0._wp 281 IF( ln_abl .AND. & 282 & ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR. & 283 & jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair ) ) THEN 284 ipka = jpka ! ABL: some fields are 3D input 285 ELSE 286 ipka = 1 287 ENDIF 288 ! 289 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,ipka) ) 290 ! 291 IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to default) 292 IF( jfpr == jp_slp ) THEN 293 sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp ! use standard pressure in Pa 294 ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN 295 sf(jfpr)%fnow(:,:,1:ipka) = 0._wp ! no precip or no snow or no surface currents 296 ELSEIF( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) THEN 297 IF( .NOT. ln_abl ) THEN 298 DEALLOCATE( sf(jfpr)%fnow ) ! deallocate as not used in this case 299 ELSE 300 sf(jfpr)%fnow(:,:,1:ipka) = 0._wp 301 ENDIF 302 ELSEIF( jfpr == jp_cc ) THEN 303 sf(jp_cc)%fnow(:,:,1:ipka) = pp_cldf 304 ELSE 305 WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr 306 CALL ctl_stop( ctmp1 ) 307 ENDIF 269 308 ELSE !-- used field --! 270 IF( ln_abl .AND. & 271 & ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR. & 272 & jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair ) ) THEN ! ABL: some fields are 3D input 273 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) ) 274 IF( sf(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) ) 275 ELSE ! others or Bulk fields are 2D fiels 276 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 277 IF( sf(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) 278 ENDIF 309 IF( sf(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) ) ! allocate array for temporal interpolation 279 310 ! 280 311 IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 ) & 281 282 312 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & 313 & ' This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' ) 283 314 ENDIF 284 315 END DO … … 327 358 WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac 328 359 WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac 329 WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac330 360 WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))' 331 361 WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12 332 362 WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15 363 WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk 364 IF(ln_crt_fbk) THEN 365 WRITE(numout,*) ' Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta' 366 WRITE(numout,*) ' Alpha rn_stau_a = ', rn_stau_a 367 WRITE(numout,*) ' Beta rn_stau_b = ', rn_stau_b 368 ENDIF 333 369 ! 334 370 WRITE(numout,*) … … 429 465 ! ! compute the surface ocean fluxes using bulk formulea 430 466 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 431 CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! <<= in 432 & sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & ! <<= in 433 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 434 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 435 & tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out 436 437 CALL blk_oce_2( sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in 438 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), & ! <<= in 439 & sf(jp_snow)%fnow(:,:,1), tsk_m, & ! <<= in 440 & zsen, zevp ) ! <=> in out 467 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), & ! <<= in 469 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 470 & sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in 471 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 472 & tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out 473 474 CALL blk_oce_2( sf(jp_tair )%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in 475 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1), & ! <<= in 476 & sf(jp_snow )%fnow(:,:,1), tsk_m, & ! <<= in 477 & zsen, zevp ) ! <=> in out 441 478 ENDIF 442 479 ! … … 470 507 471 508 472 SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi,& ! inp473 & pslp , pst , pu , pv,& ! inp474 & pqsr , pqlw ,& ! inp475 & ptsk, pssq , pcd_du, psen , pevp )! out509 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, & ! inp 510 & pslp , pst , pu , pv, & ! inp 511 & puatm, pvatm, pqsr , pqlw , & ! inp 512 & ptsk , pssq , pcd_du, psen, pevp ) ! out 476 513 !!--------------------------------------------------------------------- 477 514 !! *** ROUTINE blk_oce_1 *** … … 498 535 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 499 536 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 537 REAL(wp), INTENT(in ), DIMENSION(:,:) :: puatm ! surface current seen by the atm at T-point (i-component) [m/s] 538 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pvatm ! surface current seen by the atm at T-point (j-component) [m/s] 500 539 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr ! 501 540 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw ! … … 508 547 INTEGER :: ji, jj ! dummy loop indices 509 548 REAL(wp) :: zztmp ! local variable 549 REAL(wp) :: zstmax, zstau 550 #if defined key_cyclone 510 551 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 552 #endif 553 REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point 511 554 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 512 555 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] … … 523 566 ptsk(:,:) = pst(:,:) + rt0 ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 524 567 568 ! --- cloud cover --- ! 569 cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1) 570 525 571 ! ----------------------------------------------------------------------------- ! 526 572 ! 0 Wind components and module at T-point relative to the moving ocean ! … … 532 578 zwnd_j(:,:) = 0._wp 533 579 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 534 DO_2D_00_00 535 pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 536 pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 580 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 581 zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 582 zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 583 ! ... scalar wind at T-point (not masked) 584 wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) 585 END_2D 586 #else 587 ! ... scalar wind module at T-point (not masked) 588 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 589 wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 537 590 END_2D 538 591 #endif 539 DO_2D_00_00540 zwnd_i(ji,jj) = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) )541 zwnd_j(ji,jj) = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) )542 END_2D543 CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )544 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)545 wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) &546 & + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1)547 548 592 ! ----------------------------------------------------------------------------- ! 549 593 ! I Solar FLUX ! … … 593 637 !#LB: because AGRIF hates functions that return something else than a scalar, need to 594 638 ! use scalar version of gamma_moist() ... 595 DO_2D_11_11 596 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 597 END_2D 598 ENDIF 599 600 639 IF( ln_tpot ) THEN 640 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_zqt 642 END_2D 643 ELSE 644 ztpot = ptair(:,:) 645 ENDIF 646 ENDIF 601 647 602 648 !! Time to call the user-selected bulk parameterization for … … 627 673 628 674 END SELECT 629 675 676 IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1)) 677 IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1)) 678 IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1)) 679 !! 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=zt 681 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=zu 683 IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity " 684 IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0 685 IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu 686 630 687 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 631 688 !! ptsk and pssq have been updated!!! … … 639 696 END IF 640 697 641 !! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef.642 !! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef.643 644 IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN645 !! If zu == zt, then ensuring once for all that:646 t_zu(:,:) = ztpot(:,:)647 q_zu(:,:) = zqair(:,:)648 ENDIF649 650 651 698 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90 652 699 ! ------------------------------------------------------------- 653 700 654 701 IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 655 !! FL do we need this multiplication by tmask ... ??? 656 DO_2D_11_11 657 zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 702 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 703 zztmp = zU_zu(ji,jj) 658 704 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 659 705 pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 660 706 psen(ji,jj) = zztmp * zch_oce(ji,jj) 661 707 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 708 rhoa(ji,jj) = rho_air( ptair(ji,jj), phumi(ji,jj), pslp(ji,jj) ) 662 709 END_2D 663 710 ELSE !== BLK formulation ==! turbulent fluxes computation 664 711 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 665 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), &666 & wndm(:,:), zU_zu(:,:), pslp(:,:), &667 & taum(:,:), psen(:,:), zqla(:,:), &668 & pEvap=pevp(:,:), prhoa=rhoa(:,:) )712 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 713 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 714 & taum(:,:), psen(:,:), zqla(:,:), & 715 & pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac ) 669 716 670 717 zqla(:,:) = zqla(:,:) * tmask(:,:,1) … … 673 720 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 674 721 675 ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 676 zcd_oce = 0._wp 677 WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 678 zwnd_i = zcd_oce * zwnd_i 679 zwnd_j = zcd_oce * zwnd_j 680 681 CALL iom_put( "taum_oce", taum ) ! output wind stress module 722 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 723 IF( wndm(ji,jj) > 0._wp ) THEN 724 zztmp = taum(ji,jj) / wndm(ji,jj) 725 #if defined key_cyclone 726 ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj) 727 ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj) 728 #else 729 ztau_i(ji,jj) = zztmp * pwndi(ji,jj) 730 ztau_j(ji,jj) = zztmp * pwndj(ji,jj) 731 #endif 732 ELSE 733 ztau_i(ji,jj) = 0._wp 734 ztau_j(ji,jj) = 0._wp 735 ENDIF 736 END_2D 737 738 IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715) 739 zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp ) ! set the max value of Stau corresponding to a wind of 3 m/s (<0) 740 DO_2D( 0, 1, 0, 1 ) ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop 741 zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) ! stau (<0) must be smaller than zstmax 742 ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) ) 743 ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) ) 744 taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) ) 745 END_2D 746 ENDIF 682 747 683 748 ! ... utau, vtau at U- and V_points, resp. 684 749 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 685 ! Note th e use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves686 DO_2D _10_10687 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( z wnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) &688 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))689 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( z wnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) &690 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))750 ! Note that coastal wind stress is not used in the code... so this extra care has no effect 751 DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T 752 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj ) ) & 753 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 754 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji ,jj+1) ) & 755 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 691 756 END_2D 692 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 757 758 IF( ln_crt_fbk ) THEN 759 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. ) 760 ELSE 761 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 762 ENDIF 763 764 CALL iom_put( "taum_oce", taum ) ! output wind stress module 693 765 694 766 IF(sn_cfctl%l_prtctl) THEN … … 766 838 767 839 ! use scalar version of L_vap() for AGRIF compatibility 768 DO_2D _11_11840 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 769 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 770 842 END_2D … … 861 933 ! 862 934 INTEGER :: ji, jj ! dummy loop indices 863 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point864 935 REAL(wp) :: zootm_su ! sea-ice surface mean temperature 865 936 REAL(wp) :: zztmp1, zztmp2 ! temporary arrays … … 872 943 ! ------------------------------------------------------------ ! 873 944 ! C-grid ice dynamics : U & V-points (same as ocean) 874 DO_2D_00_00 875 zwndi_t = ( pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj ) + puice(ji,jj) ) ) 876 zwndj_t = ( pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji ,jj-1) + pvice(ji,jj) ) ) 877 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 945 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 946 wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 878 947 END_2D 879 CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. )880 948 ! 881 949 ! Make ice-atm. drag dependent on ice concentration … … 888 956 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 889 957 ENDIF 890 891 !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef. 892 !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef. 893 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 894 963 ! local scalars ( place there for vector optimisation purposes) 895 !IF (ln_abl) rhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI)896 964 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 897 965 898 966 IF( ln_blk ) THEN 899 ! ------------------------------------------------------------ ! 900 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 901 ! ------------------------------------------------------------ ! 902 ! C-grid ice dynamics : U & V-points (same as ocean) 903 DO_2D_00_00 904 putaui(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) & 905 & + rhoa(ji ,jj) * zcd_dui(ji ,jj) ) & 906 & * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 907 pvtaui(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) & 908 & + rhoa(ji,jj ) * zcd_dui(ji,jj ) ) & 909 & * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 967 ! ---------------------------------------------------- ! 968 ! Wind stress relative to nonmoving ice ( U10m ) ! 969 ! ---------------------------------------------------- ! 970 ! 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) 910 974 END_2D 911 CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 975 ! 976 DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean). 977 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 978 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 979 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) 980 putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) ) 981 pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) ) 982 END_2D 983 CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp ) 912 984 ! 913 985 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' & 914 986 & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) 915 ELSE 987 ELSE ! ln_abl 916 988 zztmp1 = 11637800.0_wp 917 989 zztmp2 = -5897.8_wp 918 DO_2D _11_11990 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 919 991 pcd_dui(ji,jj) = zcd_dui (ji,jj) 920 992 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) … … 957 1029 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 958 1030 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 959 REAL(wp) :: zfr1, zfr2 ! local variables960 1031 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature 961 1032 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice … … 966 1037 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB 967 1038 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2 1039 REAL(wp), DIMENSION(jpi,jpj) :: ztri 968 1040 !!--------------------------------------------------------------------- 969 1041 ! … … 1046 1118 evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation 1047 1119 devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT 1048 zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean1120 zevap (:,:) = emp(:,:) + tprecip(:,:) ! evaporation over ocean !LB: removed rn_efac here, correct??? 1049 1121 1050 1122 ! --- evaporation minus precipitation --- ! 1051 1123 zsnw(:,:) = 0._wp 1052 CALL ice_ thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing1124 CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing 1053 1125 emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 1054 1126 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw … … 1077 1149 END DO 1078 1150 1079 ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- ! 1080 zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! transmission when hi>10cm 1081 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 1082 ! 1083 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1084 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 1085 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 1086 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 1087 ELSEWHERE ! zero when hs>0 1088 qtr_ice_top(:,:,:) = 0._wp 1089 END WHERE 1090 ! 1091 1151 ! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- ! 1152 IF( nn_qtrice == 0 ) THEN 1153 ! formulation derived from Grenfell and Maykut (1977), where transmission rate 1154 ! 1) depends on cloudiness 1155 ! 2) is 0 when there is any snow 1156 ! 3) tends to 1 for thin ice 1157 ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:) ! surface transmission when hi>10cm 1158 DO jl = 1, jpl 1159 WHERE ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1160 qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 1161 ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp ) ! constant (ztri) when hi>10cm 1162 qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:) 1163 ELSEWHERE ! zero when hs>0 1164 qtr_ice_top(:,:,jl) = 0._wp 1165 END WHERE 1166 ENDDO 1167 ELSEIF( nn_qtrice == 1 ) THEN 1168 ! formulation is derived from the thesis of M. Lebrun (2019). 1169 ! It represents the best fit using several sets of observations 1170 ! It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90) 1171 qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:) 1172 ENDIF 1173 ! 1092 1174 IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 1093 1175 ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) … … 1171 1253 ! 1172 1254 DO jl = 1, jpl 1173 DO_2D _11_111255 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1174 1256 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness 1175 1257 IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor … … 1186 1268 ! 1187 1269 DO jl = 1, jpl 1188 DO_2D _11_111270 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1189 1271 ! 1190 1272 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness … … 1334 1416 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg] 1335 1417 ! 1336 DO_2D _00_001418 DO_2D( 0, 0, 0, 0 ) 1337 1419 ! Virtual potential temperature [K] 1338 1420 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean … … 1377 1459 ! 1378 1460 END_2D 1379 CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1. , pch, 'T', 1.)1461 CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1.0_wp, pch, 'T', 1.0_wp ) 1380 1462 ! 1381 1463 END SUBROUTINE Cdn10_Lupkes2015
Note: See TracChangeset
for help on using the changeset viewer.
