- Timestamp:
- 2014-05-15T15:56:53+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r4624 r4644 38 38 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 39 39 USE prtctl ! Print control 40 USE sbcwave,ONLY : cdn_wave !wave module40 !USE sbcwave,ONLY : cdn_wave !wave module ! moved to sbc_oce 41 41 #if defined key_lim3 || defined key_cice 42 42 USE sbc_ice ! Surface boundary condition: ice fields … … 64 64 65 65 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read) 66 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Cd_n10 ! 10m neutral drag coefficient for output 66 67 67 68 ! !!! CORE bulk parameters … … 83 84 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 84 85 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 86 LOGICAL :: ln_cdec = .FALSE. ! logical flag for using ec neutral wind drag coef 87 ! References : P. Janssen, 2008, ECMWF Workshop on Ocean-Atmosphere Interactions, 10-12 November p47-60 85 88 86 89 !! * Substitutions … … 137 140 TYPE(FLD_N) :: sn_tdif ! " " 138 141 NAMELIST/namsbc_core/ cn_dir , ln_2m , ln_taudif, rn_pfac, rn_efac, rn_vfac, & 139 & sn_wndi, sn_wndj, sn_humi , sn_qsr , &140 & sn_qlw , sn_tair, sn_prec , sn_snow, &142 & sn_wndi, sn_wndj, sn_humi , sn_qsr , ln_cdec, & 143 & sn_qlw , sn_tair, sn_prec , sn_snow, & 141 144 & sn_tdif, rn_zqt , ln_bulk2z, rn_zu 142 145 !!--------------------------------------------------------------------- … … 163 166 sn_qsr%ln_tint = .false. 164 167 ENDIF 168 IF( ln_cdec .AND. lwp ) WRITE(numout,*)'Using Hans Hersbach formula for drag.' 165 169 ! ! store namelist information in an array 166 170 slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj … … 182 186 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 183 187 ! 188 ! Drag coefficent so we can write to disk 189 ALLOCATE(Cd_n10(jpi,jpj)) 190 184 191 sfx(:,:) = 0._wp ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 185 192 ! … … 239 246 INTEGER :: ji, jj ! dummy loop indices 240 247 REAL(wp) :: zcoef_qsatw, zztmp ! local variable 248 REAL(wp) :: ztheta ! local variable, wind direction 241 249 REAL(wp), DIMENSION(:,:), POINTER :: zwnd_i, zwnd_j ! wind speed components at T-point 242 250 REAL(wp), DIMENSION(:,:), POINTER :: zqsatw ! specific humidity at pst … … 283 291 !CDIR COLLAPSE 284 292 #endif 285 DO jj = 2, jpjm1 286 DO ji = fs_2, fs_jpim1 ! vect. opt. 287 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 288 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 293 ! Wind vector at T points 294 IF ( ln_cdgw ) THEN 295 ! Use neutral 10-m wind speed from wave model if available 296 DO jj = 2, jpjm1 297 DO ji = fs_2, fs_jpim1 ! vect. opt. 298 ! Local wind direction from non-neutral 10-m wind vector (relative to grid and no current) 299 ! since we do not have the directional information from the wave model 300 ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1), sf(jp_wndj)%fnow(ji,jj,1)) 301 ! Wind vector magnitude from 10-m neutral wind speed from wave model 302 zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta) 303 zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta) 304 ! Correct for surface current, 0.0 <= rn_vfac <= 1.0 305 zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_vfac*(pu(ji-1,jj ) + pu(ji,jj))) 306 zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_vfac*(pv(ji ,jj-1) + pv(ji,jj))) 307 END DO 289 308 END DO 290 END DO 309 ELSE 310 DO jj = 2, jpjm1 311 DO ji = fs_2, fs_jpim1 ! vect. opt. 312 ! Correct for surface current, 0.0 <= rn_vfac <= 1.0 313 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - 0.5*rn_vfac*( pu(ji-1,jj ) + pu(ji,jj) ) ) 314 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - 0.5*rn_vfac*( pv(ji ,jj-1) + pv(ji,jj) ) ) 315 END DO 316 END DO 317 ENDIF 291 318 CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. ) 292 319 CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) … … 353 380 354 381 ! ... tau module, i and j component 355 DO jj = 1, jpj 356 DO ji = 1, jpi 357 zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj) 358 taum (ji,jj) = zztmp * wndm (ji,jj) 359 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 360 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 361 END DO 362 END DO 382 ! Use cdn_wave instead of Cd if we have wave parameters 383 IF ( ln_cdgw ) THEN 384 DO jj = 1, jpj 385 DO ji = 1, jpi 386 !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj) 387 zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj) 388 !taum(ji,jj) = tauoc_wavepar(ji,jj) 389 !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj) 390 taum (ji,jj) = zztmp * wndm (ji,jj) 391 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 392 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 393 ENDDO 394 ENDDO 395 ELSE 396 DO jj = 1, jpj 397 DO ji = 1, jpi 398 zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj) 399 taum (ji,jj) = zztmp * wndm (ji,jj) 400 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 401 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 402 ENDDO 403 ENDDO 404 ENDIF 363 405 364 406 ! ... add the HF tau contribution to the wind stress module? … … 379 421 CALL lbc_lnk( utau(:,:), 'U', -1. ) 380 422 CALL lbc_lnk( vtau(:,:), 'V', -1. ) 423 424 IF (ln_tauoc) THEN 425 utau(:,:) = utau(:,:)*tauoc_wavepar(:,:) 426 vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:) 427 taum(:,:) = taum(:,:)*tauoc_wavepar(:,:) 428 ENDIF 381 429 382 430 ! Turbulent fluxes over ocean … … 755 803 REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K] 756 804 REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K] 757 REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient758 805 REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient 759 806 REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient … … 775 822 ! 776 823 CALL wrk_alloc( jpi,jpj, stab ) ! integer 777 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )824 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 778 825 CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 779 826 780 827 !! * Start 781 828 !! Air/sea differences 782 dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s829 dU10 = MAX(0.5, dU) ! we don't want to fall under 0.5 m/s 783 830 dT = T_a - sst ! assuming that T_a is allready the potential temp. at zzu 784 831 dq = q_a - q_sat … … 788 835 !! 789 836 !! Neutral Drag Coefficient 790 stab = 0.5 + sign(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 791 IF ( ln_cdgw ) THEN 792 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 793 Cd_n10(:,:) = cdn_wave 837 stab = 0.5 + SIGN(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 838 839 ! Use drag coefficient from a wave model ... 840 IF (ln_cdgw) THEN 841 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 842 Cd_n10(:,:) = cdn_wave(:,:) 794 843 ELSE 795 Cd_n10 = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 844 ! ... or choose a drag law 845 IF (ln_cdec) THEN 846 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21) ! PJ (6) 847 ELSE 848 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 849 ENDIF 796 850 ENDIF 797 851 sqrt_Cd_n10 = sqrt(Cd_n10) … … 824 878 825 879 !! Updating the neutral 10m transfer coefficients : 826 Cd_n10 = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 880 IF (ln_cdec) THEN 881 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * U_n10**1.48 ) / (U_n10**0.21) ! PJ (6) 882 ELSE 883 Cd_n10 = 1E-3 * ( 2.7/U_n10 + 0.142 + U_n10/13.09 ) ! L & Y eq. (6a) 884 ENDIF 827 885 sqrt_Cd_n10 = sqrt(Cd_n10) 828 886 Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) … … 847 905 !! 848 906 CALL wrk_dealloc( jpi,jpj, stab ) ! integer 849 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )907 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 850 908 CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 851 909 ! … … 891 949 REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K] 892 950 REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K] 893 REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient894 951 REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient 895 952 REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient … … 911 968 IF( nn_timing == 1 ) CALL timing_start('TURB_CORE_2Z') 912 969 ! 913 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )970 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 914 971 CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 915 972 CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) … … 923 980 !! Neutral Drag Coefficient : 924 981 stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE 925 IF( ln_cdgw ) THEN 926 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 927 Cd_n10(:,:) = cdn_wave 982 983 ! Use drag coefficient from a wave model ... 984 IF (ln_cdgw) THEN 985 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 986 Cd_n10(:,:) = cdn_wave 928 987 ELSE 929 Cd_n10 = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 988 ! ... or choose a drag law 989 IF (ln_cdec) THEN 990 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21) ! PJ (6) 991 ELSE 992 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 993 ENDIF 930 994 ENDIF 931 995 sqrt_Cd_n10 = sqrt(Cd_n10) … … 996 1060 END DO 997 1061 !! 998 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )1062 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 999 1063 CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 1000 1064 CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
Note: See TracChangeset
for help on using the changeset viewer.
