- Timestamp:
- 2020-07-01T15:01:22+02:00 (4 years ago)
- Location:
- branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC
- Files:
-
- 3 edited
-
albedo.F90 (modified) (18 diffs)
-
sbcblk_core.F90 (modified) (34 diffs)
-
sbcssr.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
r11442 r13191 21 21 USE lib_mpp ! MPP library 22 22 USE wrk_nemo ! work arrays 23 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 23 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 24 24 USE stopack 25 25 … … 31 31 32 32 INTEGER :: albd_init = 0 !: control flag for initialization 33 33 34 34 REAL(wp) :: rmue = 0.40 ! cosine of local solar altitude 35 35 REAL(wp) :: ralb_oce = 0.066 ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) … … 37 37 REAL(wp) :: c2 = 0.10 ! " " 38 38 REAL(wp) :: rcloud = 0.06 ! cloud effect on albedo (only-for nn_ice_alb=0) 39 39 40 40 ! !!* namelist namsbc_alb 41 41 INTEGER :: nn_ice_alb … … 52 52 !!---------------------------------------------------------------------- 53 53 !! *** ROUTINE albedo_ice *** 54 !! 55 !! ** Purpose : Computation of the albedo of the snow/ice system 56 !! 54 !! 55 !! ** Purpose : Computation of the albedo of the snow/ice system 56 !! 57 57 !! ** Method : Two schemes are available (from namelist parameter nn_ice_alb) 58 58 !! 0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies … … 73 73 !! ** Note : The parameterization from Shine & Henderson-Sellers presents several misconstructions: 74 74 !! 1) ice albedo when ice thick. tends to 0 is different than ocean albedo 75 !! 2) for small ice thick. covered with some snow (<3cm?), albedo is larger 75 !! 2) for small ice thick. covered with some snow (<3cm?), albedo is larger 76 76 !! under melting conditions than under freezing conditions 77 !! 3) the evolution of ice albedo as a function of ice thickness shows 77 !! 3) the evolution of ice albedo as a function of ice thickness shows 78 78 !! 3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic 79 79 !! 80 80 !! References : Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250. 81 81 !! Brandt et al. 2005, J. Climate, vol 18 82 !! Grenfell & Perovich 2004, JGR, vol 109 82 !! Grenfell & Perovich 2004, JGR, vol 109 83 83 !!---------------------------------------------------------------------- 84 84 REAL(wp), INTENT(in ), DIMENSION(:,:,:) :: pt_ice ! ice surface temperature (Kelvin) … … 97 97 98 98 ijpl = SIZE( pt_ice, 3 ) ! number of ice categories 99 99 100 100 CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it ) 101 101 102 IF( albd_init == 0 ) CALL albedo_init ! initialization 103 104 102 IF( albd_init == 0 ) CALL albedo_init ! initialization 103 104 105 105 SELECT CASE ( nn_ice_alb ) 106 106 … … 109 109 !------------------------------------------ 110 110 CASE( 0 ) 111 111 112 112 ralb_sf = 0.80 ! dry snow 113 113 ralb_sm = 0.65 ! melting snow 114 114 ralb_if = 0.72 ! bare frozen ice 115 ralb_im = rn_albice ! bare puddled ice 116 115 ralb_im = rn_albice ! bare puddled ice 116 117 117 ! Computation of ice albedo (free of snow) 118 118 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice ) ; zalb(:,:,:) = ralb_im 119 119 ELSE WHERE ; zalb(:,:,:) = ralb_if 120 120 END WHERE 121 121 122 122 WHERE ( 1.5 < ph_ice ) ; zalb_it = zalb 123 123 ELSE WHERE( 1.0 < ph_ice .AND. ph_ice <= 1.5 ) ; zalb_it = 0.472 + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) … … 127 127 ELSE WHERE ; zalb_it = 0.1 + 3.6 * ph_ice 128 128 END WHERE 129 129 130 130 DO jl = 1, ijpl 131 131 DO jj = 1, jpj … … 133 133 ! freezing snow 134 134 ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 135 ! ! freezing snow 135 ! ! freezing snow 136 136 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 137 137 zalb_sf = ( 1._wp - zswitch ) * ( zalb_it(ji,jj,jl) & 138 138 & + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1 ) & 139 & + zswitch * ralb_sf 139 & + zswitch * ralb_sf 140 140 141 141 ! melting snow … … 143 143 zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 144 144 zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 ) & 145 & + zswitch * ralb_sm 145 & + zswitch * ralb_sm 146 146 ! 147 147 ! snow albedo 148 zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 148 zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 149 149 zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 150 150 151 151 ! Ice/snow albedo 152 152 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) … … 155 155 END DO 156 156 END DO 157 157 158 #if defined key_traldf_c2d || key_traldf_c3d 158 159 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 159 160 & CALL spp_gen( 1, pa_ice_cs(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 160 161 #else 162 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 163 & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// & 164 'key_traldf_c2d or key_traldf_c3d') 165 #endif 161 166 END DO 162 167 … … 166 171 ! New parameterization (2016) 167 172 !------------------------------------------ 168 CASE( 1 ) 173 CASE( 1 ) 169 174 170 175 ralb_im = rn_albice ! bare puddled ice … … 181 186 ! ralb_sm = 0.82 ! melting snow 182 187 ! ralb_if = 0.54 ! bare frozen ice 183 ! 188 ! 184 189 ! Computation of ice albedo (free of snow) 185 z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 190 z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 186 191 z1_c2 = 1. / 0.05 187 192 WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice ) ; zalb = ralb_im 188 193 ELSE WHERE ; zalb = ralb_if 189 194 END WHERE 190 195 191 196 WHERE ( 1.5 < ph_ice ) ; zalb_it = zalb 192 197 ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 ) ; zalb_it = zalb + ( 0.18 - zalb ) * z1_c1 * & … … 205 210 206 211 ! snow albedo 207 zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 212 zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) ) 208 213 zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 209 214 210 ! Ice/snow albedo 215 ! Ice/snow albedo 211 216 zswitch = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 212 217 pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch * zalb_it(ji,jj,jl) … … 214 219 END DO 215 220 END DO 216 221 222 #if defined key_traldf_c2d || key_traldf_c3d 217 223 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 218 224 & CALL spp_gen( 1, pa_ice_os(:,:,jl), nn_spp_icealb, rn_icealb_sd, jk_spp_alb, jl ) 219 225 #else 226 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 227 & CALL ctl_stop( 'albedo_ice: parameter perturbation will only work with '// & 228 'key_traldf_c2d or key_traldf_c3d') 229 #endif 220 230 END DO 221 231 ! Effect of the clouds (2d order polynomial) 222 pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 232 pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 223 233 224 234 END SELECT 225 235 226 236 CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it ) 227 237 ! … … 232 242 !!---------------------------------------------------------------------- 233 243 !! *** ROUTINE albedo_oce *** 234 !! 244 !! 235 245 !! ** Purpose : Computation of the albedo of the ocean 236 246 !!---------------------------------------------------------------------- … … 238 248 REAL(wp), DIMENSION(:,:), INTENT(out) :: pa_oce_cs ! albedo of ocean under clear sky 239 249 !! 240 REAL(wp) :: zcoef 250 REAL(wp) :: zcoef 241 251 !!---------------------------------------------------------------------- 242 252 ! 243 253 zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 ) ! Parameterization of Briegled and Ramanathan, 1982 244 pa_oce_cs(:,:) = zcoef 254 pa_oce_cs(:,:) = zcoef 245 255 pa_oce_os(:,:) = 0.06 ! Parameterization of Kondratyev, 1969 and Payne, 1972 246 256 ! … … 257 267 !!---------------------------------------------------------------------- 258 268 INTEGER :: ios ! Local integer output status for namelist read 259 NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 269 NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 260 270 !!---------------------------------------------------------------------- 261 271 ! -
branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r11442 r13191 161 161 ! ! ====================== ! 162 162 ! 163 rn_sfac = 1._wp ! Default to one if missing from namelist 163 rn_sfac = 1._wp ! Default to one if missing from namelist 164 164 REWIND( numnam_ref ) ! Namelist namsbc_core in reference namelist : CORE bulk parameters 165 165 READ ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901) … … 205 205 ENDIF 206 206 207 #if defined key_traldf_c2d || key_traldf_c3d 207 208 IF( ln_stopack .AND. nn_spp_relw > 0 ) THEN 208 209 rn_vfac0(:,:) = rn_vfac 209 210 CALL spp_gen(kt, rn_vfac0, nn_spp_relw, rn_relw_sd, jk_spp_relw ) 210 211 ENDIF 212 #else 213 IF ( ln_stopack .AND. nn_spp_relw > 0 ) & 214 & CALL ctl_stop( 'sbc_blk_core: parameter perturbation will only work with '// & 215 'key_traldf_c2d or key_traldf_c3d') 216 #endif 211 217 212 218 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step … … 217 223 #if defined key_cice 218 224 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 219 qlw_ice(:,:,1) = sf(jp_qlw)%fnow(:,:,1) 225 qlw_ice(:,:,1) = sf(jp_qlw)%fnow(:,:,1) 220 226 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 221 227 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 222 228 ENDIF 223 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 229 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 224 230 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 225 231 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac … … 231 237 ! 232 238 END SUBROUTINE sbc_blk_core 233 234 239 240 235 241 SUBROUTINE blk_oce_core( kt, sf, pst, pu, pv ) 236 242 !!--------------------------------------------------------------------- … … 242 248 !! ** Method : CORE bulk formulea for the ocean using atmospheric 243 249 !! fields read in sbc_read 244 !! 250 !! 245 251 !! ** Outputs : - utau : i-component of the stress at U-point (N/m2) 246 252 !! - vtau : j-component of the stress at V-point (N/m2) … … 280 286 ! local scalars ( place there for vector optimisation purposes) 281 287 zcoef_qsatw = 0.98 * 640380. / rhoa 282 288 283 289 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 284 290 … … 288 294 289 295 ! ... components ( U10m - U_oce ) at T-point (unmasked) 290 zwnd_i(:,:) = 0.e0 296 zwnd_i(:,:) = 0.e0 291 297 zwnd_j(:,:) = 0.e0 292 298 #if defined key_cyclone … … 332 338 CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm, & 333 339 & Cd, Ch, Ce, zt_zu, zq_zu ) 334 340 335 341 ! ... tau module, i and j component 336 342 DO jj = 1, jpj … … 363 369 CALL lbc_lnk( vtau(:,:), 'V', -1. ) 364 370 365 371 366 372 ! Turbulent fluxes over ocean 367 373 ! ----------------------------- … … 388 394 CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce_core: zst : ') 389 395 ENDIF 390 396 391 397 ! ----------------------------------------------------------------------------- ! 392 398 ! III Total FLUXES ! … … 396 402 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1) 397 403 ! 398 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 404 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 399 405 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip 400 406 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST … … 437 443 ! 438 444 END SUBROUTINE blk_oce_core 439 440 445 446 441 447 #if defined key_lim2 || defined key_lim3 442 448 SUBROUTINE blk_ice_core_tau … … 531 537 532 538 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_tau') 533 539 534 540 END SUBROUTINE blk_ice_core_tau 535 541 … … 544 550 !! between atmosphere and sea-ice using CORE bulk 545 551 !! formulea, ice variables and read atmmospheric fields. 546 !! 552 !! 547 553 !! caution : the net upward water flux has with mm/day unit 548 554 !!--------------------------------------------------------------------- … … 564 570 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core_flx') 565 571 ! 566 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 572 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 567 573 568 574 ! local scalars ( place there for vector optimisation purposes) … … 587 593 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 588 594 ! lw sensitivity 589 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 595 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 590 596 591 597 ! ----------------------------! … … 597 603 z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 598 604 ! Latent Heat 599 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * wndm_ice(ji,jj) & 605 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * wndm_ice(ji,jj) & 600 606 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1) ) ) 601 607 ! Latent heat sensitivity for ice (Dqla/Dt) … … 628 634 629 635 #if defined key_lim3 630 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 636 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 631 637 632 638 ! --- evaporation --- ! … … 638 644 ! --- evaporation minus precipitation --- ! 639 645 zsnw(:,:) = 0._wp 640 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing 646 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing 641 647 emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 642 648 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw … … 661 667 DO jl = 1, jpl 662 668 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 663 ! But we do not have Tice => consider it at 0 degC => evap=0 669 ! But we do not have Tice => consider it at 0 degC => evap=0 664 670 END DO 665 671 666 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 672 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 667 673 #endif 668 674 … … 688 694 ! 689 695 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_flx') 690 696 691 697 END SUBROUTINE blk_ice_core_flx 692 698 #endif … … 701 707 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 702 708 !! 703 !! ** Method : Monin Obukhov Similarity Theory 709 !! ** Method : Monin Obukhov Similarity Theory 704 710 !! + Large & Yeager (2004,2008) closure: CD_n10 = f(U_n10) 705 711 !! … … 711 717 !! - better first guess of stability by checking air-sea difference of virtual temperature 712 718 !! rather than temperature difference only... 713 !! - added function "cd_neutral_10m" that uses the improved parametrization of 719 !! - added function "cd_neutral_10m" that uses the improved parametrization of 714 720 !! Large & Yeager 2008. Drag-coefficient reduction for Cyclone conditions! 715 721 !! - using code-wide physical constants defined into "phycst.mod" rather than redifining them … … 746 752 747 753 IF( nn_timing == 1 ) CALL timing_start('turb_core_2z') 748 754 749 755 CALL wrk_alloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 750 756 CALL wrk_alloc( jpi,jpj, zeta_u, stab ) … … 758 764 U_zu = MAX( 0.5 , dU ) ! relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 759 765 760 !! First guess of stability: 766 !! First guess of stability: 761 767 ztmp0 = T_zt*(1. + 0.608*q_zt) - sst*(1. + 0.608*q_sat) ! air-sea difference of virtual pot. temp. at zt 762 768 stab = 0.5 + sign(0.5,ztmp0) ! stab = 1 if dTv > 0 => STABLE, 0 if unstable … … 772 778 Ce_n10 = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 773 779 Ch_n10 = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 774 780 775 781 !! Initializing transf. coeff. with their first guess neutral equivalents : 776 782 Cd = ztmp0 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt_Cd_n10 … … 783 789 ! 784 790 ztmp1 = T_zu - sst ! Updating air/sea differences 785 ztmp2 = q_zu - q_sat 791 ztmp2 = q_zu - q_sat 786 792 787 793 ! Updating turbulent scales : (L&Y 2004 eq. (7)) 788 794 ztmp1 = Ch/sqrt_Cd*ztmp1 ! theta* 789 795 ztmp2 = Ce/sqrt_Cd*ztmp2 ! q* 790 796 791 797 ztmp0 = T_zu*(1. + 0.608*q_zu) ! virtual potential temperature at zu 792 798 793 799 ! Estimate the inverse of Monin-Obukov length (1/L) at height zu: 794 ztmp0 = (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 800 ztmp0 = (vkarmn*grav/ztmp0*(ztmp1*(1.+0.608*q_zu) + 0.608*T_zu*ztmp2)) / (Cd*U_zu*U_zu) 795 801 ! ( Cd*U_zu*U_zu is U*^2 at zu) 796 802 … … 799 805 zpsi_h_u = psi_h( zeta_u ) 800 806 zpsi_m_u = psi_m( zeta_u ) 801 807 802 808 !! Shifting temperature and humidity at zu (L&Y 2004 eq. (9b-9c)) 803 809 IF ( .NOT. l_zt_equal_zu ) THEN … … 808 814 q_zu = max(0., q_zu) 809 815 END IF 810 816 811 817 IF( ln_cdgw ) THEN ! surface wave case 812 sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 818 sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 813 819 Cd = sqrt_Cd * sqrt_Cd 814 820 ELSE … … 820 826 ztmp0 = cd_neutral_10m(ztmp0) ! Cd_n10 821 827 sqrt_Cd_n10 = sqrt(ztmp0) 822 828 823 829 Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L&Y 2004 eq. (6b) 824 830 stab = 0.5 + sign(0.5,zeta_u) ! update stability … … 827 833 !! Update of transfer coefficients: 828 834 ztmp1 = 1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - zpsi_m_u) ! L&Y 2004 eq. (10a) 829 Cd = ztmp0 / ( ztmp1*ztmp1 ) 835 Cd = ztmp0 / ( ztmp1*ztmp1 ) 830 836 sqrt_Cd = SQRT( Cd ) 831 837 ENDIF … … 833 839 ztmp0 = (LOG(zu/10.) - zpsi_h_u) / vkarmn / sqrt_Cd_n10 834 840 ztmp2 = sqrt_Cd / sqrt_Cd_n10 835 ztmp1 = 1. + Ch_n10*ztmp0 841 ztmp1 = 1. + Ch_n10*ztmp0 836 842 Ch = Ch_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10b) 837 843 ! 838 ztmp1 = 1. + Ce_n10*ztmp0 844 ztmp1 = 1. + Ce_n10*ztmp0 839 845 Ce = Ce_n10*ztmp2 / ztmp1 ! L&Y 2004 eq. (10c) 840 846 ! … … 854 860 FUNCTION cd_neutral_10m( zw10 ) 855 861 !!---------------------------------------------------------------------- 856 !! Estimate of the neutral drag coefficient at 10m as a function 862 !! Estimate of the neutral drag coefficient at 10m as a function 857 863 !! of neutral wind speed at 10m 858 864 !! … … 860 866 !! 861 867 !! Author: L. Brodeau, june 2014 862 !!---------------------------------------------------------------------- 868 !!---------------------------------------------------------------------- 863 869 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zw10 ! scalar wind speed at 10m (m/s) 864 870 REAL(wp), DIMENSION(jpi,jpj) :: cd_neutral_10m 865 871 ! 866 872 REAL(wp), DIMENSION(:,:), POINTER :: rgt33 867 !!---------------------------------------------------------------------- 873 !!---------------------------------------------------------------------- 868 874 ! 869 875 CALL wrk_alloc( jpi,jpj, rgt33 ) 870 876 ! 871 877 !! When wind speed > 33 m/s => Cyclone conditions => special treatment 872 rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) ) ! If zw10 < 33. => 0, else => 1 878 rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) ) ! If zw10 < 33. => 0, else => 1 873 879 cd_neutral_10m = 1.e-3 * ( & 874 880 & (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. -
branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
r11442 r13191 24 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 25 USE timing ! Timing 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 27 USE stopack 28 28 USE wrk_nemo ! Memory Allocation … … 42 42 REAL(wp) :: rn_dqdt ! restoring factor on SST and SSS 43 43 REAL(wp) :: rn_deds ! restoring factor on SST and SSS 44 LOGICAL :: ln_sssr_bnd ! flag to bound erp term 44 LOGICAL :: ln_sssr_bnd ! flag to bound erp term 45 45 REAL(wp) :: rn_sssr_bnd ! ABS(Max./Min.) value of erp term [mm/day] 46 46 … … 101 101 rn_dqdt_s(:,:) = rn_dqdt 102 102 103 IF( ln_stopack .AND. nn_spp_dqdt > 0 ) & 104 & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt ) 103 #if defined key_traldf_c2d || key_traldf_c3d 104 IF( ln_stopack .AND. nn_spp_dqdt > 0 ) & 105 & CALL spp_gen( kt, rn_dqdt_s, nn_spp_dqdt, rn_dqdt_sd, jk_spp_dqdt ) 106 #else 107 IF ( ln_stopack .AND. nn_spp_dqdt > 0 ) & 108 & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 109 'key_traldf_c2d or key_traldf_c3d') 110 #endif 111 105 112 DO jj = 1, jpj 106 113 DO ji = 1, jpi … … 117 124 CALL wrk_alloc( jpi, jpj, zsrp) 118 125 zsrp(:,:) = rn_deds 126 #if defined key_traldf_c2d || key_traldf_c3d 119 127 IF( ln_stopack .AND. nn_spp_dedt > 0 ) & 120 128 & CALL spp_gen(kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 129 #else 130 IF ( ln_stopack .AND. nn_spp_icealb > 0 ) & 131 & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 132 'key_traldf_c2d or key_traldf_c3d') 133 #endif 134 135 121 136 !CDIR COLLAPSE 122 137 DO jj = 1, jpj 123 138 DO ji = 1, jpi 124 139 zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) ) & ! No damping in vicinity of river mouths 125 & * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 140 & * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) 126 141 sfx(ji,jj) = sfx(ji,jj) + zerp ! salt flux 127 142 erp(ji,jj) = zerp / MAX( sss_m(ji,jj), 1.e-20 ) ! converted into an equivalent volume flux (diagnostic only) … … 134 149 CALL wrk_alloc( jpi, jpj, zsrp) 135 150 zsrp(:,:) = rn_deds 151 #if defined key_traldf_c2d || key_traldf_c3d 136 152 IF( ln_stopack .AND. nn_spp_dedt > 0 ) & 137 153 & CALL spp_gen( kt, zsrp, nn_spp_dedt, rn_dedt_sd, jk_spp_deds ) 138 zerp_bnd = rn_sssr_bnd / rday ! - - 154 #else 155 IF ( ln_stopack .AND. nn_spp_dedt > 0 ) & 156 & CALL ctl_stop( 'sbc_ssr: parameter perturbation will only work with '// & 157 'key_traldf_c2d or key_traldf_c3d') 158 #endif 159 zerp_bnd = rn_sssr_bnd / rday ! - - 139 160 !CDIR COLLAPSE 140 161 DO jj = 1, jpj 141 DO ji = 1, jpi 162 DO ji = 1, jpi 142 163 zerp = (zsrp(ji,jj)/rday) * ( 1. - 2.*rnfmsk(ji,jj) ) & ! No damping in vicinity of river mouths 143 164 & * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) & … … 161 182 END SUBROUTINE sbc_ssr 162 183 163 184 164 185 SUBROUTINE sbc_ssr_init 165 186 !!--------------------------------------------------------------------- … … 184 205 !!---------------------------------------------------------------------- 185 206 ! 186 187 REWIND( numnam_ref ) ! Namelist namsbc_ssr in reference namelist : 207 208 REWIND( numnam_ref ) ! Namelist namsbc_ssr in reference namelist : 188 209 READ ( numnam_ref, namsbc_ssr, IOSTAT = ios, ERR = 901) 189 210 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_ssr in reference namelist', lwp ) … … 240 261 ENDIF 241 262 ! 242 ! !* Initialize qrp and erp if no restoring 263 ! !* Initialize qrp and erp if no restoring 243 264 IF( nn_sstr /= 1 ) qrp(:,:) = 0._wp 244 265 IF( nn_sssr /= 1 .OR. nn_sssr /= 2 ) erp(:,:) = 0._wp 245 266 ! 246 267 END SUBROUTINE sbc_ssr_init 247 268 248 269 !!====================================================================== 249 270 END MODULE sbcssr
Note: See TracChangeset
for help on using the changeset viewer.
