- Timestamp:
- 2020-06-10T13:13:39+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r11442 r13088 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.
Note: See TracChangeset
for help on using the changeset viewer.