- Timestamp:
- 2020-11-27T17:26:33+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_icb_1900
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/SBC/sbcblk.F90 (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_icb_1900
- Property svn:externals
-
NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcblk.F90
r13226 r13899 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) … … 87 87 INTEGER , PUBLIC, PARAMETER :: jp_voatm = 11 ! index of surface current (j-component) 88 88 ! ! seen by the atmospheric forcing (m/s) at T-point 89 INTEGER , PUBLIC, PARAMETER :: jp_hpgi = 12 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 90 INTEGER , PUBLIC, PARAMETER :: jp_hpgj = 13 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 91 INTEGER , PUBLIC, PARAMETER :: jpfld = 13 ! maximum number of files to read 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 92 93 93 94 ! Warning: keep this structure allocatable for Agrif... … … 175 176 TYPE(FLD_N) :: sn_qlw , sn_tair , sn_prec, sn_snow ! " " 176 177 TYPE(FLD_N) :: sn_slp , sn_uoatm, sn_voatm ! " " 177 TYPE(FLD_N) :: sn_ hpgi, sn_hpgj! " "178 TYPE(FLD_N) :: sn_cc, sn_hpgi, sn_hpgj ! " " 178 179 INTEGER :: ipka ! number of levels in the atmospheric variable 179 180 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 180 181 & sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm, & 181 & sn_ hpgi, sn_hpgj,&182 & sn_cc, sn_hpgi, sn_hpgj, & 182 183 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 183 184 & cn_dir , rn_zqt, rn_zu, & … … 260 261 slf_i(jp_tair ) = sn_tair ; slf_i(jp_humi ) = sn_humi 261 262 slf_i(jp_prec ) = sn_prec ; slf_i(jp_snow ) = sn_snow 262 slf_i(jp_slp ) = sn_slp 263 slf_i(jp_slp ) = sn_slp ; slf_i(jp_cc ) = sn_cc 263 264 slf_i(jp_uoatm) = sn_uoatm ; slf_i(jp_voatm) = sn_voatm 264 265 slf_i(jp_hpgi ) = sn_hpgi ; slf_i(jp_hpgj ) = sn_hpgj … … 289 290 ! 290 291 IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to default) 291 IF( jfpr == jp_slp 292 IF( jfpr == jp_slp ) THEN 292 293 sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp ! use standard pressure in Pa 293 294 ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN 294 295 sf(jfpr)%fnow(:,:,1:ipka) = 0._wp ! no precip or no snow or no surface currents 295 ELSEIF( ( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) .AND. .NOT. ln_abl ) THEN 296 DEALLOCATE( sf(jfpr)%fnow ) ! deallocate as not used in this case 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 297 304 ELSE 298 305 WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr … … 303 310 ! 304 311 IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 ) & 305 306 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...' ) 307 314 ENDIF 308 315 END DO … … 559 566 ptsk(:,:) = pst(:,:) + rt0 ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 560 567 568 ! --- cloud cover --- ! 569 cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1) 570 561 571 ! ----------------------------------------------------------------------------- ! 562 572 ! 0 Wind components and module at T-point relative to the moving ocean ! … … 568 578 zwnd_j(:,:) = 0._wp 569 579 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 570 DO_2D _11_11580 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 571 581 zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 572 582 zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) … … 576 586 #else 577 587 ! ... scalar wind module at T-point (not masked) 578 DO_2D _11_11588 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 579 589 wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 580 590 END_2D … … 628 638 ! use scalar version of gamma_moist() ... 629 639 IF( ln_tpot ) THEN 630 DO_2D _11_11640 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 631 641 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 632 642 END_2D … … 690 700 691 701 IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 692 DO_2D _11_11702 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 693 703 zztmp = zU_zu(ji,jj) 694 704 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod … … 710 720 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 711 721 712 DO_2D _11_11722 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 713 723 IF( wndm(ji,jj) > 0._wp ) THEN 714 724 zztmp = taum(ji,jj) / wndm(ji,jj) … … 728 738 IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715) 729 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) 730 DO_2D _01_01! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop740 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 731 741 zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) ! stau (<0) must be smaller than zstmax 732 742 ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) ) … … 739 749 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 740 750 ! Note that coastal wind stress is not used in the code... so this extra care has no effect 741 DO_2D _00_00! start loop at 2, in case ln_crt_fbk = T751 DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T 742 752 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj ) ) & 743 753 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) … … 828 838 829 839 ! use scalar version of L_vap() for AGRIF compatibility 830 DO_2D _11_11840 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 831 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 832 842 END_2D … … 933 943 ! ------------------------------------------------------------ ! 934 944 ! C-grid ice dynamics : U & V-points (same as ocean) 935 DO_2D _11_11945 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 936 946 wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 937 947 END_2D … … 959 969 ! ---------------------------------------------------- ! 960 970 ! supress moving ice in wind stress computation as we don't know how to do it properly... 961 DO_2D _01_01! at T point971 DO_2D( 0, 1, 0, 1 ) ! at T point 962 972 putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndi(ji,jj) 963 973 pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndj(ji,jj) 964 974 END_2D 965 975 ! 966 DO_2D _00_00! U & V-points (same as ocean).976 DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean). 967 977 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 968 978 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) … … 978 988 zztmp1 = 11637800.0_wp 979 989 zztmp2 = -5897.8_wp 980 DO_2D _11_11990 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 981 991 pcd_dui(ji,jj) = zcd_dui (ji,jj) 982 992 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) … … 1019 1029 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 1020 1030 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 1021 REAL(wp) :: zfr1, zfr2 ! local variables1022 1031 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature 1023 1032 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice … … 1028 1037 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB 1029 1038 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2 1039 REAL(wp), DIMENSION(jpi,jpj) :: ztri 1030 1040 !!--------------------------------------------------------------------- 1031 1041 ! … … 1112 1122 ! --- evaporation minus precipitation --- ! 1113 1123 zsnw(:,:) = 0._wp 1114 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 1115 1125 emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 1116 1126 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw … … 1139 1149 END DO 1140 1150 1141 ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- ! 1142 zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! transmission when hi>10cm 1143 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 1144 ! 1145 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1146 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 1147 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 1148 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 1149 ELSEWHERE ! zero when hs>0 1150 qtr_ice_top(:,:,:) = 0._wp 1151 END WHERE 1152 ! 1153 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 ! 1154 1174 IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 1155 1175 ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) … … 1233 1253 ! 1234 1254 DO jl = 1, jpl 1235 DO_2D _11_111255 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1236 1256 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness 1237 1257 IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor … … 1248 1268 ! 1249 1269 DO jl = 1, jpl 1250 DO_2D _11_111270 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1251 1271 ! 1252 1272 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness … … 1396 1416 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg] 1397 1417 ! 1398 DO_2D _00_001418 DO_2D( 0, 0, 0, 0 ) 1399 1419 ! Virtual potential temperature [K] 1400 1420 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean
Note: See TracChangeset
for help on using the changeset viewer.
