- Timestamp:
- 2016-07-01T18:02:45+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r5602 r6772 44 44 USE sbc_ice ! Surface boundary condition: ice fields 45 45 USE lib_fortran ! to use key_nosignedzero 46 USE sbcapr 46 47 #if defined key_lim3 47 48 USE ice, ONLY : u_ice, v_ice, jpl, pfrld, a_i_b 48 49 USE limthd_dh ! for CALL lim_thd_snwblow 49 50 #elif defined key_lim2 50 USE ice_2, ONLY : u_ice, v_ice 51 USE ice_2, ONLY : u_ice, v_ice, pfrld 51 52 USE par_ice_2 52 53 #endif … … 83 84 REAL(wp), PARAMETER :: Cice = 1.4e-3 ! iovi 1.63e-3 ! transfer coefficient over ice 84 85 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 86 REAL(wp), PARAMETER :: rgas = 287.1 ! gas const. dry air (J/kg/K) 87 REAL(wp), PARAMETER :: rvap = 461.51 ! gas const. vapour (J/kg/K) 85 88 86 89 ! !!* Namelist namsbc_core : CORE bulk parameters … … 91 94 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 92 95 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 96 ! 97 LOGICAL :: ln_tair_celsius !: logical flag for Read Tair: Tair in NEMO is Kelvin 98 LOGICAL :: ln_humi_rel !: logical flag for Read relative humidity (T) or specific humidity (F) 99 LOGICAL :: ln_cohum_arc !: logical flag for Correction of Humidity in the Arctic Ocean 100 LOGICAL :: ln_cotair_arc !: logical flag for Correction of Air Temperature in the Arctic Ocean 101 LOGICAL :: ln_corad_antar !: logical flag for Correction of radiatives fluxes in the Southern Ocean 102 93 103 94 104 !! * Substitutions … … 143 153 INTEGER :: ios ! Local integer output status for namelist read 144 154 ! 155 INTEGER :: ji,jj 156 REAL(wp) :: zzlat, zzlat1, zzlat2, zfm, zfrld 157 REAL(wp) :: zmin,zmax 158 REAL(wp), DIMENSION(:,:), POINTER :: xyt,z_qsr,z_qlw,z_qsr1,z_qlw1, z_hum, z_tair 159 REAL(wp), DIMENSION(:,:), POINTER :: zqsr_lr, zqsr_hr, zqlw_lr, zqlw_hr 160 145 161 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 146 162 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read … … 151 167 & sn_wndi, sn_wndj, sn_humi , sn_qsr , & 152 168 & sn_qlw , sn_tair, sn_prec , sn_snow, & 153 & sn_tdif, rn_zqt, rn_zu 169 & sn_tdif, rn_zqt, rn_zu , ln_tair_celsius, & 170 & ln_humi_rel , ln_cohum_arc, & 171 & ln_cotair_arc, ln_corad_antar 172 154 173 !!--------------------------------------------------------------------- 155 174 ! 175 CALL wrk_alloc( jpi,jpj, xyt,z_qsr,z_qlw,z_qsr1,z_qlw1, z_hum, z_tair ) 176 CALL wrk_alloc( jpi,jpj, zqsr_lr, zqsr_hr, zqlw_lr, zqlw_hr ) 156 177 ! ! ====================== ! 157 178 IF( kt == nit000 ) THEN ! First call kt=nit000 ! … … 194 215 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 195 216 ! 217 ! 218 IF(lwp) WRITE(numout,*) 'sbc_blk_core: jfld = ',jfld 219 IF( ln_cohum_arc ) CALL ctl_warn( 'sbc_blk_core: correction of humidity in arctic' ) 220 IF( ln_cotair_arc ) CALL ctl_warn( 'sbc_blk_core: correction of air temperature in arctic' ) 221 IF( ln_corad_antar ) CALL ctl_warn( 'sbc_blk_core: correction of short radiation in antartic' ) 222 IF( ln_humi_rel ) CALL ctl_warn( 'sbc_blk_core: use relative humidity instead of specific humidity') 223 IF( ln_tair_celsius) CALL ctl_warn( 'sbc_blk_core: Tair is read in Celsius') 224 IF(lwp) WRITE(numout,*) 'sbc_blk_core: rn_pfac = ',rn_pfac 225 ! 196 226 sfx(:,:) = 0._wp ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 197 227 ! … … 199 229 200 230 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 231 232 !========================================= 233 ! ONLINE CORRECTIONS 234 !========================================= 235 ! 236 ! Correction of Tair 237 ! 238 IF( ln_tair_celsius .AND. MOD( kt-1, nn_fsbc ) == 0 ) THEN 239 sf(jp_tair)%fnow = sf(jp_tair)%fnow + 273.15_wp ! Conversion of the Temperature °C --> Kelvin 240 ENDIF 241 ! 242 ! Correction of SW and LW in the Southern Ocean 243 ! 244 IF( ln_corad_antar .AND. .NOT. sf(jp_qsr)%ln_tint .AND. MOD( kt-1, 86400/INT(rdt) ) == 0 ) THEN 245 z_qsr(:,:) = 0.8 * sf(jp_qsr)%fnow(:,:,1) 246 xyt(:,:) = 0.e0 ; zzlat1 = -65. ; zzlat2 = -60. 247 DO jj = 1, jpj 248 DO ji = 1, jpi 249 zzlat = gphit(ji,jj) 250 IF( zzlat >= zzlat1 .AND. zzlat <= zzlat2 ) THEN 251 xyt(ji,jj) = (zzlat2-zzlat)/(zzlat2-zzlat1) 252 ELSE IF ( zzlat < zzlat1 ) THEN 253 xyt(ji,jj) = 1 254 ENDIF 255 END DO 256 END DO 257 IF(lwp) WRITE(numout,*) 'Correc ln_corad_antar' 258 z_qsr1(:,:) = z_qsr(:,:) * xyt(:,:) + ( 1.0 - xyt(:,:) ) * sf(jp_qsr)%fnow(:,:,1) 259 sf(jp_qsr)%fnow(:,:,1) = z_qsr1(:,:) 260 ENDIF 261 262 IF( MOD( kt-1, nn_fsbc ) == 0 )THEN 263 ! 264 IF ( nmonth >= 5 .AND. nmonth <= 9 ) THEN 265 ! 266 ! Correction of Humidity in the Arctic Ocean 267 ! 268 IF( ln_cohum_arc ) THEN 269 z_hum(:,:) = 0.85 * sf(jp_humi)%fnow(:,:,1) 270 xyt(:,:) = 0.e0 ; zzlat1 = 78. ; zzlat2 = 82. 271 DO jj = 1, jpj 272 DO ji = 1, jpi 273 zzlat = gphit(ji,jj) 274 #if defined key_lim2 || defined key_lim3 275 IF ( ALLOCATED(pfrld) ) THEN ; zfrld = pfrld(ji,jj) ; ELSE ; zfrld = 0 ; ENDIF 276 #endif 277 IF( zzlat >= zzlat1 .AND. zzlat <= zzlat2 .AND. zfrld < 0.85 ) THEN 278 xyt(ji,jj) = ( zzlat - zzlat1 ) / ( zzlat2 - zzlat1 ) 279 ELSE IF ( zzlat > zzlat2 .AND. zfrld < 0.85 ) THEN 280 xyt(ji,jj) = 1._wp 281 ENDIF 282 ENDDO 283 ENDDO 284 IF(lwp) WRITE(numout,*) 'Correc ln_cohum_arc' 285 sf(jp_humi)%fnow(:,:,1) = z_hum(:,:) * xyt(:,:) + ( 1.0 - xyt(:,:) ) * sf(jp_humi)%fnow(:,:,1) 286 ENDIF 287 ! 288 ! Correction of Air Temperature in the Arctic Ocean 289 ! 290 IF( ln_cotair_arc ) THEN 291 z_tair(:,:) = sf(jp_tair)%fnow(:,:,1) - 2.0 292 xyt(:,:) = 0.e0 ; zzlat1 = 78. ; zzlat2 = 82. 293 DO jj = 1, jpj 294 DO ji = 1, jpi 295 zzlat = gphit(ji,jj) 296 #if defined key_lim2 || defined key_lim3 297 IF( ALLOCATED(pfrld) ) THEN ; zfrld = pfrld(ji,jj) ; ELSE ; zfrld=0 ; ENDIF 298 #endif 299 IF( zzlat >= zzlat1 .AND. zzlat <= zzlat2 .AND. zfrld < 0.85 ) THEN 300 xyt(ji,jj) = ( zzlat - zzlat1 ) / ( zzlat2 - zzlat1 ) 301 ELSE IF( zzlat > zzlat2 .AND. zfrld < 0.85 ) THEN 302 xyt(ji,jj) = 1._wp 303 ENDIF 304 END DO 305 ENDDO 306 IF(lwp) WRITE(numout,*) 'Correc ln_cotair_arc' 307 sf(jp_tair)%fnow(:,:,1) = z_tair(:,:) * xyt(:,:) + ( 1.0 - xyt(:,:) ) * sf(jp_tair)%fnow(:,:,1) 308 ENDIF 309 ! 310 ENDIF ! 5 <= nmonth <= 9 311 312 ! 313 ENDIF ! IF MOD( kt-1, nn_fsbc ) 314 315 DO jj=1,jpj 316 DO ji=1,jpi 317 sf(jp_humi)%fnow(ji,jj,1) = MAX( MIN( sf(jp_humi)%fnow(ji,jj,1) ,1.0 ) , 0.0 ) 318 sf(jp_prec)%fnow(ji,jj,1) = MAX( sf(jp_prec)%fnow(ji,jj,1) ,0.0 ) 319 sf(jp_qsr )%fnow(ji,jj,1) = MAX( sf(jp_qsr )%fnow(ji,jj,1) ,0.0 ) 320 sf(jp_qlw )%fnow(ji,jj,1) = MAX( sf(jp_qlw )%fnow(ji,jj,1) ,0.0 ) 321 ENDDO 322 END DO 323 324 ! 325 !========================================= 326 ! END OF ONLINE CORRECTIONS 327 !========================================= 328 ! 201 329 202 330 ! ! compute the surface ocean fluxes using CORE bulk formulea … … 215 343 ENDIF 216 344 #endif 345 ! 346 CALL wrk_dealloc( jpi,jpj, xyt,z_qsr,z_qlw,z_qsr1,z_qlw1, z_hum, z_tair ) 347 CALL wrk_dealloc( jpi,jpj, zqsr_lr, zqsr_hr, zqlw_lr, zqlw_hr ) 217 348 ! 218 349 END SUBROUTINE sbc_blk_core … … 257 388 REAL(wp), DIMENSION(:,:), POINTER :: zt_zu ! air temperature at wind speed height 258 389 REAL(wp), DIMENSION(:,:), POINTER :: zq_zu ! air spec. hum. at wind speed height 390 REAL(wp), DIMENSION(:,:), POINTER :: zqatm , zpatm ! specific humidity and mean sea level pressure (Pa) 391 REAL(wp) :: vt, vp, vq, zqa, zq0, zq1, zq2, zee 259 392 !!--------------------------------------------------------------------- 260 393 ! … … 262 395 ! 263 396 CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap ) 264 CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )397 CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu ,zqatm, zpatm ) 265 398 ! 266 399 ! local scalars ( place there for vector optimisation purposes) … … 314 447 ! ... specific humidity at SST and IST 315 448 zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 316 449 ! 450 IF ( ln_humi_rel ) THEN 451 zq0 = rvap / rgas - 1.0 452 zq1 = rgas / rvap 453 zq2 = 1.0 - zq1 454 zpatm(:,:) = 100800. ! atmospheric pressure (assumed constant here) 455 IF ( ln_apr_dyn ) zpatm(:,:) = apr(:,:) 456 DO jj = 1 , jpj 457 DO ji = 1 , jpi 458 vt = sf(jp_tair)%fnow(ji,jj,1) - rt0 ! air temperature (Celsius) 459 vp = zpatm(ji,jj) / 100. ! mean sea level pressure (mb or hPa) 460 vq = sf(jp_humi)%fnow(ji,jj,1) ! relative humidity (fraction of 1) 461 ! Convert RH at the air/sea interface in specific humidity (kg/kg) 462 ! Teten's formula for qsat (mb) 463 zqa = ( 1.0007 + 3.46e-6 * vp) * 6.1121 * EXP( 17.502 * vt / ( 240.97+vt ) ) 464 zee = zqa * vq ! vapour partial pressure (mb) 465 vq = zq1 * zee / ( vp - zq2 * zee ) ! specific humidity (kg/kg) 466 zqatm(ji,jj) = vq 467 ENDDO 468 ENDDO 469 ELSE 470 zqatm(:,:)=sf(jp_humi)%fnow(:,:,1) 471 ENDIF 472 ! 317 473 ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point : 318 CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm, &474 CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, zqatm, wndm, & 319 475 & Cd, Ch, Ce, zt_zu, zq_zu ) 320 476 … … 354 510 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 355 511 !! q_air and t_air are (or "are almost") given at 10m (wind reference height) 356 zevap(:,:) = rn_efac*MAX( 0._wp, rhoa*Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) )*wndm(:,:) ) ! Evaporation 512 !zevap(:,:) = rn_efac*MAX( 0._wp, rhoa*Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) )*wndm(:,:) ) ! Evaporation 513 zevap(:,:) = rn_efac*MAX( 0._wp, rhoa*Ce(:,:)*( zqsatw(:,:) - zqatm(:,:) )*wndm(:,:) ) ! Evaporation 357 514 zqsb (:,:) = cpa*rhoa*Ch(:,:)*( zst (:,:) - sf(jp_tair)%fnow(:,:,1) )*wndm(:,:) ! Sensible Heat 358 515 ELSE … … 414 571 ! 415 572 CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap ) 416 CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu )573 CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu, zqatm, zpatm ) 417 574 ! 418 575 IF( nn_timing == 1 ) CALL timing_stop('blk_oce_core') … … 437 594 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 438 595 !!--------------------------------------------------------------------- 439 ! 596 440 597 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core_tau') 441 598 ! … … 530 687 REAL(wp) :: zcoef_dqlw, zcoef_dqla, zcoef_dqsb 531 688 REAL(wp) :: zztmp, z1_lsub ! temporary variable 689 REAL(wp) :: ztamr,zmt1,zmt2,zmt3,zev,zes 532 690 !! 533 691 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice … … 536 694 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 537 695 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 696 REAL(wp), DIMENSION(:,:) , POINTER :: zqatm, zpatm , ztatm ! specific humidity 538 697 !!--------------------------------------------------------------------- 539 698 ! … … 541 700 ! 542 701 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 702 CALL wrk_alloc( jpi,jpj, zqatm, zpatm, ztatm ) 703 704 IF ( ln_humi_rel ) THEN 705 zpatm(:,:) = 100800. ! atmospheric pressure (assumed constant here) 706 IF (ln_apr_dyn) zpatm(:,:) = apr(:,:) 707 DO jj=1,jpj 708 DO ji=1,jpi 709 ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj,1) ! air temperature in Kelvins 710 ztamr = ztatm(ji,jj) - rtt ! Saturation water vapour 711 zmt1 = SIGN( 17.269, ztamr ) 712 zmt2 = SIGN( 21.875, ztamr ) 713 zmt3 = SIGN( 28.200, -ztamr ) 714 zes = 611.0 * EXP( ABS( ztamr ) * MIN ( zmt1, zmt2 ) & 715 & / ( ztatm(ji,jj) - 35.86 + MAX( 0.e0, zmt3 ) ) ) 716 zev = sf(jp_humi)%fnow(ji,jj,1) * zes ! vapour pressure 717 zqatm(ji,jj) = 0.622 * zev / ( zpatm(ji,jj) - 0.378 * zev ) ! specific humidity 718 ENDDO 719 ENDDO 720 ELSE 721 zqatm(:,:) = sf(jp_humi)%fnow(:,:,1) 722 ENDIF 543 723 544 724 ! local scalars ( place there for vector optimisation purposes) … … 574 754 ! Latent Heat 575 755 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * wndm_ice(ji,jj) & 576 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1) ) )756 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - zqatm(ji,jj) ) ) 577 757 ! Latent heat sensitivity for ice (Dqla/Dt) 578 758 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN … … 659 839 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_flx') 660 840 841 CALL wrk_dealloc( jpi,jpj, zqatm, zpatm, ztatm ) 661 842 END SUBROUTINE blk_ice_core_flx 662 843 #endif
Note: See TracChangeset
for help on using the changeset viewer.
