- Timestamp:
- 2018-11-07T18:25:49+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9866_HPC_03_globcom/src/OCE/SBC/sbcblk.F90
- Property svn:keywords set to Id
r9767 r10288 46 46 USE lib_fortran ! to use key_nosignedzero 47 47 #if defined key_si3 48 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su, rn_cnd_s 48 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su, rn_cnd_s, hfx_err_dif 49 49 USE icethd_dh ! for CALL ice_thd_snwblow 50 50 #endif … … 135 135 !!---------------------------------------------------------------------- 136 136 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 137 !! $Id : sbcblk.F90 6416 2016-04-01 12:22:17Z clem$138 !! Software governed by the CeCILL licen ce (./LICENSE)137 !! $Id$ 138 !! Software governed by the CeCILL license (see ./LICENSE) 139 139 !!---------------------------------------------------------------------- 140 140 CONTAINS … … 239 239 !drag coefficient read from wave model definable only with mfs bulk formulae and core 240 240 ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR ) THEN 241 CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')241 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 242 242 ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN 243 243 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') … … 526 526 ! 527 527 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 528 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus& ! remove latent melting heat for solid precip528 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 529 529 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST 530 530 & + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair 531 531 & * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & 532 532 & + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 533 & * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 _snow ) - rt0 ) * cpic533 & * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi 534 534 qns(:,:) = qns(:,:) * tmask(:,:,1) 535 535 ! … … 659 659 zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 660 660 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT 661 gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) )661 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 662 662 END DO 663 663 END DO … … 792 792 REAL(wp) :: zst3 ! local variable 793 793 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 794 REAL(wp) :: zztmp, z1_ lsub ! - -794 REAL(wp) :: zztmp, z1_rLsub ! - - 795 795 REAL(wp) :: zfr1, zfr2 ! local variables 796 796 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature … … 868 868 869 869 ! --- evaporation --- ! 870 z1_ lsub = 1._wp /Lsub871 evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_ lsub ! sublimation872 devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_ lsub ! d(sublimation)/dT873 zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean870 z1_rLsub = 1._wp / rLsub 871 evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation 872 devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT 873 zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean 874 874 875 875 ! --- evaporation minus precipitation --- ! … … 884 884 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair 885 885 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) 886 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 _snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )886 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 887 887 qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only) 888 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 _snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )888 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 889 889 890 890 ! --- total solar and non solar fluxes --- ! … … 894 894 895 895 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 896 qprec_ice(:,:) = rhos n * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )896 qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 897 897 898 898 ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 899 899 DO jl = 1, jpl 900 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic* tmask(:,:,1) )900 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 901 901 ! ! But we do not have Tice => consider it at 0degC => evap=0 902 902 END DO … … 907 907 ! 908 908 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 909 q sr_ice_tr(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )909 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 910 910 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 911 q sr_ice_tr(:,:,:) = qsr_ice(:,:,:) * zfr1911 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 912 912 ELSEWHERE ! zero when hs>0 913 q sr_ice_tr(:,:,:) = 0._wp913 qtr_ice_top(:,:,:) = 0._wp 914 914 END WHERE 915 915 ! … … 971 971 CASE ( 1 , 2 ) 972 972 ! 973 zfac = 1._wp / ( rn_cnd_s + rc dic)973 zfac = 1._wp / ( rn_cnd_s + rcnd_i ) 974 974 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 975 975 zfac3 = 2._wp / zepsilon … … 978 978 DO jj = 1 , jpj 979 979 DO ji = 1, jpi 980 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rc dic * phs(ji,jj,jl) ) * zfac! Effective thickness980 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness 981 981 IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 982 982 END DO … … 990 990 ! -------------------------------------------------------------! 991 991 ! 992 zfac = rc dic* rn_cnd_s992 zfac = rcnd_i * rn_cnd_s 993 993 ! 994 994 DO jl = 1, jpl … … 997 997 ! 998 998 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 999 & ( rc dic* phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )999 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 1000 1000 ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature 1001 1001 ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature 1002 zqa0 = qsr_ice(ji,jj,jl) - q sr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl)! Net initial atmospheric heat flux1002 zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 1003 1003 ! 1004 1004 DO iter = 1, nit ! --- Iterative loop … … 1011 1011 qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 1012 1012 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1013 qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - q sr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )&1013 qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) ) & 1014 1014 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1015 1016 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1017 hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 1015 1018 1016 1019 END DO
Note: See TracChangeset
for help on using the changeset viewer.
