- Timestamp:
- 2014-01-16T15:27:29+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_ECMWF_waves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r4333 r4351 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 … … 138 141 NAMELIST/namsbc_core/ cn_dir , ln_2m , ln_taudif, rn_pfac, rn_efac, rn_vfac, & 139 142 & sn_wndi, sn_wndj, sn_humi , sn_qsr , & 143 & sn_ccov, ln_cdec, ln_cldcov, rn_rrelwind 140 144 & sn_qlw , sn_tair, sn_prec , sn_snow, & 141 145 & sn_tdif, rn_zqt , ln_bulk2z, rn_zu … … 163 167 sn_qsr%ln_tint = .false. 164 168 ENDIF 169 IF( ln_cdec .AND. lwp ) WRITE(numout,*)'Using Hans Hersbach formula for drag.' 165 170 ! ! store namelist information in an array 166 171 slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj … … 182 187 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 183 188 ! 189 ! Drag coefficent so we can write to disk 190 ALLOCATE(Cd_n10(jpi,jpj)) 191 184 192 sfx(:,:) = 0._wp ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 185 193 ! … … 239 247 INTEGER :: ji, jj ! dummy loop indices 240 248 REAL(wp) :: zcoef_qsatw, zztmp ! local variable 249 REAL(wp) :: ztheta ! local variable, wind direction 241 250 REAL(wp), DIMENSION(:,:), POINTER :: zwnd_i, zwnd_j ! wind speed components at T-point 242 251 REAL(wp), DIMENSION(:,:), POINTER :: zqsatw ! specific humidity at pst … … 283 292 !CDIR COLLAPSE 284 293 #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) ) ) 294 ! Wind vector at T points 295 IF ( ln_cdgw ) THEN 296 ! Use neutral 10-m wind speed from wave model if available 297 DO jj = 2, jpjm1 298 DO ji = fs_2, fs_jpim1 ! vect. opt. 299 ! Local wind direction from non-neutral 10-m wind vector (relative to grid and no current) 300 ! since we do not have the directional information from the wave model 301 ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1), sf(jp_wndj)%fnow(ji,jj,1)) 302 ! Wind vector magnitude from 10-m neutral wind speed from wave model 303 zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta) 304 zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta) 305 ! Correct for surface current, 0.0 <= rn_rrelwind <= 1.0 306 zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_rrelwind*(pu(ji-1,jj ) + pu(ji,jj))) 307 zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_rrelwind*(pv(ji ,jj-1) + pv(ji,jj))) 308 END DO 289 309 END DO 290 END DO 310 ELSE 311 DO jj = 2, jpjm1 312 DO ji = fs_2, fs_jpim1 ! vect. opt. 313 ! Correct for surface current, 0.0 <= rn_rrelwind <= 1.0 314 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - 0.5*rn_rrelwind*( pu(ji-1,jj ) + pu(ji,jj) ) ) 315 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - 0.5*rn_rrelwind*( pv(ji ,jj-1) + pv(ji,jj) ) ) 316 END DO 317 END DO 318 ENDIF 291 319 CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. ) 292 320 CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) … … 353 381 354 382 ! ... 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 383 ! Use cdn_wave instead of Cd if we have wave parameters 384 IF ( ln_cdgw ) THEN 385 DO jj = 1, jpj 386 DO ji = 1, jpi 387 !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj) 388 zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj) 389 !taum(ji,jj) = tauoc_wavepar(ji,jj) 390 !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj) 391 taum (ji,jj) = zztmp * wndm (ji,jj) 392 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 393 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 394 ENDDO 395 ENDDO 396 ELSE 397 DO jj = 1, jpj 398 DO ji = 1, jpi 399 zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj) 400 taum (ji,jj) = zztmp * wndm (ji,jj) 401 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 402 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 403 ENDDO 404 ENDDO 405 ENDIF 363 406 364 407 ! ... add the HF tau contribution to the wind stress module? … … 379 422 CALL lbc_lnk( utau(:,:), 'U', -1. ) 380 423 CALL lbc_lnk( vtau(:,:), 'V', -1. ) 424 425 IF (ln_tauoc) THEN 426 utau(:,:) = utau(:,:)*tauoc_wavepar(:,:) 427 vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:) 428 taum(:,:) = taum(:,:)*tauoc_wavepar(:,:) 429 ENDIF 381 430 382 431 ! Turbulent fluxes over ocean … … 755 804 REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K] 756 805 REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K] 757 REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient758 806 REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient 759 807 REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient … … 775 823 ! 776 824 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 )825 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 778 826 CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 779 827 780 828 !! * Start 781 829 !! Air/sea differences 782 dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s830 dU10 = MAX(0.5, dU) ! we don't want to fall under 0.5 m/s 783 831 dT = T_a - sst ! assuming that T_a is allready the potential temp. at zzu 784 832 dq = q_a - q_sat … … 788 836 !! 789 837 !! 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 838 stab = 0.5 + SIGN(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 839 840 ! Use drag coefficient from a wave model ... 841 IF (ln_cdgw) THEN 842 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 843 Cd_n10(:,:) = cdn_wave(:,:) 794 844 ELSE 795 Cd_n10 = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 845 ! ... or choose a drag law 846 IF (ln_cdec) THEN 847 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21) ! PJ (6) 848 ELSE 849 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 850 ENDIF 796 851 ENDIF 797 852 sqrt_Cd_n10 = sqrt(Cd_n10) … … 824 879 825 880 !! 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) 881 IF (ln_cdec) THEN 882 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * U_n10**1.48 ) / (U_n10**0.21) ! PJ (6) 883 ELSE 884 Cd_n10 = 1E-3 * ( 2.7/U_n10 + 0.142 + U_n10/13.09 ) ! L & Y eq. (6a) 885 ENDIF 827 886 sqrt_Cd_n10 = sqrt(Cd_n10) 828 887 Ce_n10 = 1.e-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) … … 847 906 !! 848 907 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 )908 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 850 909 CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 851 910 ! … … 891 950 REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K] 892 951 REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K] 893 REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient894 952 REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient 895 953 REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient … … 911 969 IF( nn_timing == 1 ) CALL timing_start('TURB_CORE_2Z') 912 970 ! 913 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )971 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 914 972 CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 915 973 CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) … … 923 981 !! Neutral Drag Coefficient : 924 982 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 983 984 ! Use drag coefficient from a wave model ... 985 IF (ln_cdgw) THEN 986 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 987 Cd_n10(:,:) = cdn_wave 928 988 ELSE 929 Cd_n10 = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 989 ! ... or choose a drag law 990 IF (ln_cdec) THEN 991 Cd_n10 = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21) ! PJ (6) 992 ELSE 993 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 994 ENDIF 930 995 ENDIF 931 996 sqrt_Cd_n10 = sqrt(Cd_n10) … … 996 1061 END DO 997 1062 !! 998 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, C d_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )1063 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 999 1064 CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 1000 1065 CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
Note: See TracChangeset
for help on using the changeset viewer.
