Changeset 12098 for NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE
- Timestamp:
- 2019-12-06T17:12:55+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90
r12074 r12098 24 24 !! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition 25 25 !! sbc_blk : bulk formulation as ocean surface boundary condition 26 !! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE) 27 !! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE) 26 !! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE) 27 !! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE) 28 28 !! sea-ice case only : 29 29 !! blk_ice_1 : provide the air-ice stress … … 209 209 !! Some namelist sanity tests: 210 210 IF( ln_NCAR ) & 211 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' ) 211 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' ) 212 212 IF( nn_fsbc /= 1 ) & 213 & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.') 213 & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.') 214 214 END IF 215 215 216 216 IF( ln_skin_wl ) THEN 217 217 !! Check if the frequency of downwelling solar flux input makes sense and if ln_dm2dc=T if it is daily! … … 221 221 & CALL ctl_stop( 'sbc_blk_init: Please set ln_dm2dc=T for warm-layer param. (ln_skin_wl) to work properly' ) 222 222 END IF 223 223 224 224 ioptio = 0 225 225 IF( ln_humi_sph ) THEN … … 254 254 slf_i(jp_slp ) = sn_slp 255 255 IF( ln_abl ) THEN 256 slf_i(jp_hpgi) = sn_hpgi ; slf_i(jp_hpgj) = sn_hpgj256 slf_i(jp_hpgi) = sn_hpgi ; slf_i(jp_hpgj) = sn_hpgj 257 257 END IF 258 258 ! … … 263 263 DO jfpr= 1, jpfld 264 264 ! 265 IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to zero) 265 IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to zero) 266 266 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 267 sf(jfpr)%fnow(:,:,1) = 0._wp 267 sf(jfpr)%fnow(:,:,1) = 0._wp 268 268 ELSE !-- used field --! 269 269 IF( ln_abl .AND. & … … 304 304 ENDIF 305 305 ! 306 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient 306 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient 307 307 rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level 308 308 rn_zu = ght_abl(2) … … 364 364 !! (momentum, heat, freshwater and runoff) 365 365 !! 366 !! ** Method : 366 !! ** Method : 367 367 !! (1) READ each fluxes in NetCDF files: 368 368 !! the wind velocity (i-component) at z=rn_zu (m/s) at T-point … … 374 374 !! the total precipitation (rain+snow) (Kg/m2/s) 375 375 !! the snow (solid precipitation) (kg/m2/s) 376 !! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) 376 !! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) 377 377 !! (2) CALL blk_oce_1 and blk_oce_2 378 378 !! … … 402 402 IF( kt == nit000 ) THEN 403 403 WRITE(numout,*) '' 404 #if defined key_agrif 405 WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 406 #else 404 407 ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain 405 408 IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points! … … 422 425 END IF 423 426 WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 427 #endif 424 428 WRITE(numout,*) '' 425 429 END IF !IF( kt == nit000 ) 426 427 430 ! ! compute the surface ocean fluxes using bulk formulea 428 431 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN … … 432 435 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 433 436 & zssq, zcd_du, zsen, zevp ) ! =>> out 434 437 435 438 CALL blk_oce_2( sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in 436 439 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), & ! <<= in … … 469 472 470 473 SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, & ! inp 471 472 473 474 & pslp , pst , pu , pv, & ! inp 475 & pqsr , pqlw , & ! inp 476 & pssq , pcd_du, psen , pevp ) ! out 474 477 !!--------------------------------------------------------------------- 475 478 !! *** ROUTINE blk_oce_1 *** … … 496 499 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 497 500 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 498 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr ! 499 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw ! 501 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr ! 502 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw ! 500 503 REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] 501 504 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] … … 604 607 !! Time to call the user-selected bulk parameterization for 605 608 !! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more... 606 SELECT CASE( nblk ) 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 609 SELECT CASE( nblk ) 610 611 CASE( np_NCAR ) 612 CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, & 613 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 614 615 CASE( np_COARE_3p0 ) 616 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 617 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 618 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 619 620 CASE( np_COARE_3p6 ) 621 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 622 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 623 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 624 625 CASE( np_ECMWF ) 626 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 627 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 628 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 629 630 CASE DEFAULT 628 631 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 629 632 … … 645 648 !! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef. 646 649 !! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef. 647 650 648 651 IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 649 652 !! If zu == zt, then ensuring once for all that: … … 651 654 q_zu(:,:) = zqair(:,:) 652 655 ENDIF 653 656 654 657 655 658 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90 656 659 ! ------------------------------------------------------------- 657 660 658 661 IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 659 !! FL do we need this multiplication by tmask ... ???662 !! FL do we need this multiplication by tmask ... ??? 660 663 DO jj = 1, jpj 661 664 DO ji = 1, jpi 662 665 zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 663 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 666 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 664 667 pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 665 668 psen(ji,jj) = zztmp * zch_oce(ji,jj) 666 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 669 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 667 670 END DO 668 671 END DO … … 678 681 taum(:,:) = taum(:,:) * tmask(:,:,1) 679 682 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 680 683 681 684 ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 682 685 zcd_oce = 0._wp 683 686 WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 684 687 zwnd_i = zcd_oce * zwnd_i 685 zwnd_j = zcd_oce * zwnd_j 688 zwnd_j = zcd_oce * zwnd_j 686 689 687 690 CALL iom_put( "taum_oce", taum ) ! output wind stress module 688 691 689 692 ! ... utau, vtau at U- and V_points, resp. 690 693 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines … … 718 721 719 722 SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec, & ! <<= in 720 723 & psnow, pst , psen, pevp ) ! <<= in 721 724 !!--------------------------------------------------------------------- 722 725 !! *** ROUTINE blk_oce_2 *** 723 726 !! 724 !! ** Purpose : finalize the momentum, heat and freshwater fluxes computation 725 !! at the ocean surface at each time step knowing Cd, Ch, Ce and 727 !! ** Purpose : finalize the momentum, heat and freshwater fluxes computation 728 !! at the ocean surface at each time step knowing Cd, Ch, Ce and 726 729 !! atmospheric variables (from ABL or external data) 727 730 !! … … 829 832 END SUBROUTINE blk_oce_2 830 833 831 834 832 835 #if defined key_si3 833 836 !!---------------------------------------------------------------------- … … 856 859 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s] 857 860 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s] 858 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: phumi ! atmospheric wind at T-point [m/s] 861 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: phumi ! atmospheric wind at T-point [m/s] 859 862 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s] 860 863 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! " … … 865 868 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl 866 869 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl 867 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl 870 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl 868 871 ! 869 872 INTEGER :: ji, jj ! dummy loop indices … … 894 897 Ce_ice(:,:) = Cd_ice(:,:) 895 898 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 896 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 899 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 897 900 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 898 901 ENDIF … … 904 907 !IF (ln_abl) rhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI) 905 908 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 906 907 IF( ln_blk ) THEN 909 910 IF( ln_blk ) THEN 908 911 ! ------------------------------------------------------------ ! 909 912 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! … … 926 929 ELSE 927 930 zztmp1 = 11637800.0_wp 928 931 zztmp2 = -5897.8_wp 929 932 DO jj = 1, jpj 930 933 DO ji = 1, jpi … … 1082 1085 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 1083 1086 qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 1084 1087 1085 1088 ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 1086 1089 DO jl = 1, jpl … … 1256 1259 ! ice-atm drag 1257 1260 pcd(:,:) = rCd_ice + & ! pure ice drag 1258 1261 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1259 1262 1260 1263 END SUBROUTINE Cdn10_Lupkes2012 … … 1365 1368 ! Momentum Transfert Coefficients (Eq. 38) 1366 1369 pcd(ji,jj) = zCdn_skin_ice * zfmi + & 1367 1370 & zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 1368 1371 1369 1372 ! Heat Transfert Coefficients (Eq. 49) 1370 1373 pch(ji,jj) = zChn_skin_ice * zfhi + & 1371 1374 & zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 1372 1375 ! 1373 1376 END DO
Note: See TracChangeset
for help on using the changeset viewer.
