- Timestamp:
- 2020-09-15T12:49:18+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/src/OCE/SBC/sbcblk.F90
r13467 r13469 408 408 #if defined key_cyclone 409 409 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 410 DO jj = 2, jpjm1 411 DO ji = fs_2, fs_jpim1 ! vect. opt. 412 sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj) 413 sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj) 414 END DO 415 END DO 410 DO_2D_00_00 411 sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj) 412 sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj) 413 END_2D 416 414 #endif 417 DO jj = 2, jpjm1 418 DO ji = fs_2, fs_jpim1 ! vect. opt. 419 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 420 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 421 END DO 422 END DO 415 DO_2D_00_00 416 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 417 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 418 END_2D 423 419 CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 424 420 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) … … 474 470 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 475 471 476 DO jj = 1, jpj ! tau module, i and j component 477 DO ji = 1, jpi 478 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 479 taum (ji,jj) = zztmp * wndm (ji,jj) 480 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 481 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 482 END DO 483 END DO 472 DO_2D_11_11 473 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 474 taum (ji,jj) = zztmp * wndm (ji,jj) 475 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 476 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 477 END_2D 484 478 485 479 ! ! add the HF tau contribution to the wind stress module … … 491 485 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 492 486 ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 493 DO jj = 1, jpjm1 494 DO ji = 1, fs_jpim1 495 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & 496 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 497 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) & 498 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 499 END DO 500 END DO 487 DO_2D_10_10 488 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & 489 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 490 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) & 491 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 492 END_2D 501 493 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 502 494 … … 633 625 !!---------------------------------------------------------------------------------- 634 626 ! 635 DO jj = 1, jpj 636 DO ji = 1, jpi 627 DO_2D_11_11 628 ! 629 ztmp = rt0 / ptak(ji,jj) 630 ! 631 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 632 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) & 633 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) & 634 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 ) 637 635 ! 638 ztmp = rt0 / ptak(ji,jj) 639 ! 640 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 641 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) & 642 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) & 643 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 ) 644 ! 645 q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat ) ! 0.01 because SLP is in [Pa] 646 ! 647 END DO 648 END DO 636 q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat ) ! 0.01 because SLP is in [Pa] 637 ! 638 END_2D 649 639 ! 650 640 END FUNCTION q_sat … … 669 659 !!---------------------------------------------------------------------------------- 670 660 ! 671 DO jj = 1, jpj 672 DO ji = 1, jpi 673 zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 674 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT 675 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 676 END DO 677 END DO 661 DO_2D_11_11 662 zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 663 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT 664 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 665 END_2D 678 666 ! 679 667 END FUNCTION gamma_moist … … 735 723 ! ------------------------------------------------------------ ! 736 724 ! C-grid ice dynamics : U & V-points (same as ocean) 737 DO jj = 2, jpjm1 738 DO ji = fs_2, fs_jpim1 ! vect. opt. 739 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 740 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 741 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 742 END DO 743 END DO 725 DO_2D_00_00 726 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 727 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 728 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 729 END_2D 744 730 CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. ) 745 731 ! … … 763 749 ! ------------------------------------------------------------ ! 764 750 zztmp1 = rn_vfac * 0.5_wp 765 DO jj = 2, jpj ! at T point 766 DO ji = 2, jpi 767 zztmp2 = zrhoa(ji,jj) * Cd_atm(ji,jj) * wndm_ice(ji,jj) 768 utau_ice(ji,jj) = zztmp2 * ( sf(jp_wndi)%fnow(ji,jj,1) - zztmp1 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 769 vtau_ice(ji,jj) = zztmp2 * ( sf(jp_wndj)%fnow(ji,jj,1) - zztmp1 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 770 END DO 771 END DO 772 ! 773 DO jj = 2, jpjm1 ! U & V-points (same as ocean). 774 DO ji = fs_2, fs_jpim1 ! vect. opt. 775 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 776 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 777 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) 778 utau_ice(ji,jj) = zztmp1 * ( utau_ice(ji,jj) + utau_ice(ji+1,jj ) ) 779 vtau_ice(ji,jj) = zztmp2 * ( vtau_ice(ji,jj) + vtau_ice(ji ,jj+1) ) 780 END DO 781 END DO 751 DO_2D_01_01 752 zztmp2 = zrhoa(ji,jj) * Cd_atm(ji,jj) * wndm_ice(ji,jj) 753 utau_ice(ji,jj) = zztmp2 * ( sf(jp_wndi)%fnow(ji,jj,1) - zztmp1 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 754 vtau_ice(ji,jj) = zztmp2 * ( sf(jp_wndj)%fnow(ji,jj,1) - zztmp1 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 755 END_2D 756 ! 757 DO_2D_00_00 758 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 759 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 760 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) 761 utau_ice(ji,jj) = zztmp1 * ( utau_ice(ji,jj) + utau_ice(ji+1,jj ) ) 762 vtau_ice(ji,jj) = zztmp2 * ( vtau_ice(ji,jj) + vtau_ice(ji ,jj+1) ) 763 END_2D 782 764 CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. ) 783 765 ! … … 1025 1007 ! 1026 1008 DO jl = 1, jpl 1027 DO jj = 1 , jpj 1028 DO ji = 1, jpi 1029 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness 1030 IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 1031 END DO 1032 END DO 1009 DO_2D_11_11 1010 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness 1011 IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 1012 END_2D 1033 1013 END DO 1034 1014 ! … … 1042 1022 ! 1043 1023 DO jl = 1, jpl 1044 DO jj = 1 , jpj 1045 DO ji = 1, jpi 1046 ! 1047 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1048 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 1049 ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature 1050 ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature 1051 zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 1052 ! 1053 DO iter = 1, nit ! --- Iterative loop 1054 zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards) 1055 zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget 1056 ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update 1057 END DO 1058 ! 1059 ptsu (ji,jj,jl) = MIN( rt0, ztsu ) 1060 qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 1061 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1062 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) ) & 1063 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1064 1065 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1066 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) 1067 1024 DO_2D_11_11 1025 ! 1026 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1027 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 1028 ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature 1029 ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature 1030 zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 1031 ! 1032 DO iter = 1, nit ! --- Iterative loop 1033 zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards) 1034 zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget 1035 ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update 1068 1036 END DO 1069 END DO 1037 ! 1038 ptsu (ji,jj,jl) = MIN( rt0, ztsu ) 1039 qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 1040 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1041 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) ) & 1042 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1043 1044 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1045 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) 1046 1047 END_2D 1070 1048 ! 1071 1049 END DO … … 1195 1173 zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg] 1196 1174 ! 1197 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility 1198 DO ji = fs_2, fs_jpim1 1199 ! Virtual potential temperature [K] 1200 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1201 zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1202 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1203 1204 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1205 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1206 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1207 1208 ! Momentum and Heat Neutral Transfert Coefficients 1209 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1210 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1211 1212 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 1213 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1214 z0i = z0_skin_ice ! over ice (cf Lupkes email for details) 1215 IF( zrib_o <= 0._wp ) THEN 1216 zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 10 1217 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26 1218 & )**zgamma )**z1_gamma 1219 ELSE 1220 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12 1221 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1222 ENDIF 1223 1224 IF( zrib_i <= 0._wp ) THEN 1225 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 1226 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25 1227 ELSE 1228 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11 1229 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1230 ENDIF 1231 1232 ! Momentum Transfert Coefficients (Eq. 38) 1233 Cd(ji,jj) = zCdn_skin_ice * zfmi + & 1234 & 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) ) 1235 1236 ! Heat Transfert Coefficients (Eq. 49) 1237 Ch(ji,jj) = zChn_skin_ice * zfhi + & 1238 & 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) ) 1239 ! 1240 END DO 1241 END DO 1175 DO_2D_00_00 1176 ! Virtual potential temperature [K] 1177 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1178 zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1179 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1180 1181 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1182 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1183 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1184 1185 ! Momentum and Heat Neutral Transfert Coefficients 1186 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1187 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1188 1189 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 1190 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1191 z0i = z0_skin_ice ! over ice (cf Lupkes email for details) 1192 IF( zrib_o <= 0._wp ) THEN 1193 zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 10 1194 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26 1195 & )**zgamma )**z1_gamma 1196 ELSE 1197 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12 1198 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1199 ENDIF 1200 1201 IF( zrib_i <= 0._wp ) THEN 1202 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 1203 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25 1204 ELSE 1205 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11 1206 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1207 ENDIF 1208 1209 ! Momentum Transfert Coefficients (Eq. 38) 1210 Cd(ji,jj) = zCdn_skin_ice * zfmi + & 1211 & 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) ) 1212 1213 ! Heat Transfert Coefficients (Eq. 49) 1214 Ch(ji,jj) = zChn_skin_ice * zfhi + & 1215 & 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) ) 1216 ! 1217 END_2D 1242 1218 CALL lbc_lnk_multi( 'sbcblk', Cd, 'T', 1., Ch, 'T', 1. ) 1243 1219 !
Note: See TracChangeset
for help on using the changeset viewer.
