! $Id: surf_ocean_mod.F90 5144 2024-07-29 21:01:04Z evignon $ MODULE surf_ocean_mod IMPLICIT NONE CONTAINS !****************************************************************************** SUBROUTINE surf_ocean(rlon, rlat, swnet, lwnet, alb1, & windsp, rmu0, fder, tsurf_in, & itime, dtime, jour, knon, knindex, & p1lay, z1lay, cdragh, cdragm, precip_rain, precip_snow, precip_bs, temp_air, spechum, & AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ps, u1, v1, gustiness, rugoro, pctsrf, & snow, qsurf, agesno, & z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, & tsurf_new, dflux_s, dflux_l, lmt_bils, & flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, & dt_ds, tkt, tks, taur, sss & #ifdef ISO ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, & xtsnow,xtevap,h1 & #endif ) USE albedo, ONLY: alboc, alboc_cd USE bulk_flux_m, ONLY: bulk_flux USE dimphy, ONLY: klon, zmasq USE surface_data, ONLY: type_ocean USE ocean_forced_mod, ONLY: ocean_forced_noice USE ocean_slab_mod, ONLY: ocean_slab_noice USE ocean_cpl_mod, ONLY: ocean_cpl_noice USE indice_sol_mod, ONLY: nbsrf, is_oce USE lmdz_abort_physic, ONLY: abort_physic #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso,niso #ifdef ISOVERIF USE isotopes_mod, ONLY: iso_eau,ridicule USE isotopes_verif_mod #endif #endif USE limit_read_mod USE config_ocean_skin_m, ONLY: activate_ocean_skin USE lmdz_clesphys, ONLY: iflag_cycle_diurne, iflag_z0_oce, nsw, f_z0qh_oce, fmagic, ok_bs, & iflag_albedo, pmagic, z0min ! for cycle_diurne and for iflag_z0_oce==-1 (prescribed z0) USE lmdz_yomcst ! This SUBROUTINE will make a CALL to ocean_XXX_noice according to the ocean mode (force, ! slab or couple). The calculations of albedo and rugosity for the ocean surface are ! done in here because they are identical for the different modes of ocean. ! Input variables !****************************************************************************** INTEGER, INTENT(IN) :: itime, jour, knon INTEGER, DIMENSION(klon), INTENT(IN) :: knindex REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat REAL, DIMENSION(klon), INTENT(IN) :: swnet ! net shortwave radiation at surface REAL, DIMENSION(klon), INTENT(IN) :: lwnet ! net longwave radiation at surface REAL, DIMENSION(klon), INTENT(IN) :: alb1 ! albedo in visible SW interval REAL, DIMENSION(klon), INTENT(IN) :: windsp ! wind at 10 m, in m s-1 REAL, DIMENSION(klon), INTENT(IN) :: rmu0 REAL, DIMENSION(klon), INTENT(IN) :: fder REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! defined only for subscripts 1:knon REAL, DIMENSION(klon), INTENT(IN) :: p1lay, z1lay ! pression (Pa) et altitude (m) du premier niveau REAL, DIMENSION(klon), INTENT(IN) :: cdragh REAL, DIMENSION(klon), INTENT(IN) :: cdragm REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow, precip_bs REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV REAL, DIMENSION(klon), INTENT(IN) :: ps REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness REAL, DIMENSION(klon), INTENT(IN) :: rugoro REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf #ifdef ISO REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow REAL, DIMENSION(ntraciso,klon), INTENT(IN) :: xtspechum #endif ! In/Output variables !****************************************************************************** REAL, DIMENSION(klon), INTENT(INOUT) :: snow REAL, DIMENSION(klon), INTENT(INOUT) :: qsurf REAL, DIMENSION(klon), INTENT(INOUT) :: agesno REAL, DIMENSION(klon), INTENT(inOUT) :: z0h #ifdef ISO REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsnow REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce #endif REAL, INTENT(INOUT) :: delta_sst(:) ! (knon) ! Ocean-air interface temperature minus bulk SST, in K. Defined ! only if activate_ocean_skin >= 1. REAL, INTENT(INOUT) :: delta_sal(:) ! (knon) ! Ocean-air interface salinity minus bulk salinity, in ppt. Defined ! only if activate_ocean_skin >= 1. REAL, INTENT(INOUT) :: ds_ns(:) ! (knon) ! "delta salinity near surface". Salinity variation in the ! near-surface turbulent layer. That is subskin salinity minus ! foundation salinity. In ppt. REAL, INTENT(INOUT) :: dt_ns(:) ! (knon) ! "delta temperature near surface". Temperature variation in the ! near-surface turbulent layer. That is subskin temperature ! minus foundation temperature. (Can be negative.) In K. REAL, INTENT(INOUT) :: dter(:) ! (knon) ! Temperature variation in the diffusive microlayer, that is ! ocean-air interface temperature minus subskin temperature. In ! K. REAL, INTENT(INOUT) :: dser(:) ! (knon) ! Salinity variation in the diffusive microlayer, that is ! ocean-air interface salinity minus subskin salinity. In ppt. REAL, INTENT(INOUT) :: dt_ds(:) ! (knon) ! (tks / tkt) * dTer, in K ! Output variables !************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: z0m !albedo SB >>> ! REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval ! REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval REAL, DIMENSION(6), INTENT(IN) :: SFRWL REAL, DIMENSION(klon, nsw), INTENT(OUT) :: alb_dir_new, alb_dif_new !albedo SB <<< REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! sea surface temperature, in K REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 REAL, INTENT(OUT) :: tkt(:) ! (knon) ! épaisseur (m) de la couche de diffusion thermique (microlayer) ! cool skin thickness REAL, INTENT(OUT) :: tks(:) ! (knon) ! épaisseur (m) de la couche de diffusion de masse (microlayer) REAL, INTENT(OUT) :: taur(:) ! (knon) ! momentum flux due to rain, in Pa REAL, INTENT(OUT) :: sss(:) ! (klon) ! Bulk salinity of the surface layer of the ocean, in ppt. (Only ! defined for subscripts 1:knon, but we have to declare it with ! size klon because of the coupling machinery.) #ifdef ISO REAL, DIMENSION(ntraciso,klon), INTENT(OUT) :: xtevap ! isotopes in surface evaporation flux REAL, DIMENSION(klon), INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation #endif ! Local variables !************************************************************************* INTEGER :: i, k REAL :: tmp REAL, PARAMETER :: cepdu2 = (0.1)**2 REAL, DIMENSION(klon) :: alb_eau, z0_lim REAL, DIMENSION(klon) :: radsol REAL, DIMENSION(klon) :: cdragq ! Cdrag pour l'evaporation REAL, DIMENSION(klon) :: precip_totsnow CHARACTER(len = 20), PARAMETER :: modname = "surf_ocean" REAL rhoa(knon) ! density of moist air (kg / m3) REAL sens_prec_liq(knon) REAL t_int(knon) ! ocean-air interface temperature, in K REAL s_int(knon) ! ocean-air interface salinity, in ppt !************************************************************************** #ifdef ISO #ifdef ISOVERIF DO i = 1, knon IF (iso_eau > 0) THEN CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), & spechum(i),'surf_ocean_mod 117', & errmax,errmaxrel) CALL iso_verif_egalite_choix(xtsnow(iso_eau,i), & snow(i),'surf_ocean_mod 127', & errmax,errmaxrel) ENDIF !IF (iso_eau > 0) THEN ENDDO !DO i=1,klon #endif #endif !****************************************************************************** ! Calculate total net radiance at surface !****************************************************************************** radsol(1:klon) = 0.0 ! initialisation a priori inutile radsol(1:knon) = swnet(1:knon) + lwnet(1:knon) !**************************************************************************************** !Total solid precip IF (ok_bs) THEN precip_totsnow(:) = precip_snow(:) + precip_bs(:) ELSE precip_totsnow(:) = precip_snow(:) ENDIF !****************************************************************************** ! Cdragq computed from cdrag ! The difference comes only from a factor (f_z0qh_oce) on z0, so that ! it can be computed inside surf_ocean ! More complicated appraches may require the propagation through ! pbl_surface of an independant cdragq variable. !****************************************************************************** IF (f_z0qh_oce /= 1.) THEN ! Si on suit les formulations par exemple de Tessel, on ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55 cdragq(1:knon) = cdragh(1:knon) * & log(z1lay(1:knon) / z0h(1:knon)) / log(z1lay(1:knon) / (f_z0qh_oce * z0h(1:knon))) ELSE cdragq(1:knon) = cdragh(1:knon) ENDIF rhoa = PS(:KNON) / (Rd * temp_air(:knon) * (1. + retv * spechum(:knon))) !****************************************************************************** ! Switch according to type of ocean (couple, slab or forced) !****************************************************************************** SELECT CASE(type_ocean) CASE('couple') CALL ocean_cpl_noice(& swnet, lwnet, alb1, & windsp, fder, & itime, dtime, knon, knindex, & p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, temp_air, spechum, & AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ps, u1, v1, gustiness, tsurf_in, & radsol, snow, agesno, & qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & tsurf_new, dflux_s, dflux_l, sens_prec_liq, sss, delta_sal, rhoa, & delta_sst, dTer, dSer, dt_ds) CASE('slab') CALL ocean_slab_noice(& itime, dtime, jour, knon, knindex, & p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, temp_air, spechum, & AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ps, u1, v1, gustiness, tsurf_in, & radsol, snow, & qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & tsurf_new, dflux_s, dflux_l, lmt_bils) CASE('force') CALL ocean_forced_noice(& itime, dtime, jour, knon, knindex, & p1lay, cdragh, cdragq, cdragm, precip_rain, precip_totsnow, & temp_air, spechum, & AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ps, u1, v1, gustiness, tsurf_in, & radsol, snow, agesno, & qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa & #ifdef ISO ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, & xtsnow,xtevap,h1 & #endif ) END SELECT !****************************************************************************** ! fcodron: compute lmt_bils forced case (same as wfbils_oce / 1.-contfracatm) !****************************************************************************** IF (type_ocean/='slab') THEN lmt_bils(1:klon) = 0. DO i = 1, knon lmt_bils(knindex(i)) = (swnet(i) + lwnet(i) + fluxsens(i) + fluxlat(i)) & * pctsrf(knindex(i), is_oce) / (1. - zmasq(knindex(i))) END DO END IF !****************************************************************************** ! Calculate ocean surface albedo !****************************************************************************** !albedo SB >>> IF (iflag_albedo==0) THEN !--old parametrizations of ocean surface albedo IF (iflag_cycle_diurne>=1) THEN CALL alboc_cd(rmu0, alb_eau) !--ad-hoc correction for model radiative balance tuning !--now outside alboc_cd routine alb_eau(1:klon) = fmagic * alb_eau(1:klon) + pmagic alb_eau(1:klon) = MIN(MAX(alb_eau(1:klon), 0.0), 1.0) ELSE CALL alboc(REAL(jour), rlat, alb_eau) !--ad-hoc correction for model radiative balance tuning !--now outside alboc routine alb_eau(1:klon) = fmagic * alb_eau(1:klon) + pmagic alb_eau(1:klon) = MIN(MAX(alb_eau(1:klon), 0.04), 0.60) ENDIF DO i = 1, knon DO k = 1, nsw alb_dir_new(i, k) = alb_eau(knindex(i)) ENDDO ENDDO !IM 09122015 next line corresponds to the old way of doing in LMDZ5A/IPSLCM5A versions !albedo for diffuse radiation is taken the same as for direct radiation alb_dif_new(1:knon, :) = alb_dir_new(1:knon, :) !IM 09122015 end ELSE IF (iflag_albedo==1) THEN !--new parametrization of ocean surface albedo by Sunghye Baek !--albedo for direct and diffuse radiation are different CALL ocean_albedo(knon, rmu0, knindex, windsp, SFRWL, alb_dir_new, alb_dif_new) !--ad-hoc correction for model radiative balance tuning alb_dir_new(1:knon, :) = fmagic * alb_dir_new(1:knon, :) + pmagic alb_dif_new(1:knon, :) = fmagic * alb_dif_new(1:knon, :) + pmagic alb_dir_new(1:knon, :) = MIN(MAX(alb_dir_new(1:knon, :), 0.0), 1.0) alb_dif_new(1:knon, :) = MIN(MAX(alb_dif_new(1:knon, :), 0.0), 1.0) ELSE IF (iflag_albedo==2) THEN ! F. Codron albedo read from limit.nc CALL limit_read_rug_alb(itime, dtime, jour, & knon, knindex, z0_lim, alb_eau) DO i = 1, knon DO k = 1, nsw alb_dir_new(i, k) = alb_eau(i) ENDDO ENDDO alb_dif_new = alb_dir_new ENDIF !albedo SB <<< !****************************************************************************** ! Calculate the rugosity !****************************************************************************** IF (iflag_z0_oce==0) THEN DO i = 1, knon tmp = MAX(cepdu2, gustiness(i) + u1(i)**2 + v1(i)**2) z0m(i) = 0.018 * cdragm(i) * (gustiness(i) + u1(i)**2 + v1(i)**2) / RG & + 0.11 * 14e-6 / SQRT(cdragm(i) * tmp) z0m(i) = MAX(1.5e-05, z0m(i)) ENDDO z0h(1:knon) = z0m(1:knon) ! En attendant mieux ELSE IF (iflag_z0_oce==1) THEN DO i = 1, knon tmp = MAX(cepdu2, gustiness(i) + u1(i)**2 + v1(i)**2) z0m(i) = 0.018 * cdragm(i) * (gustiness(i) + u1(i)**2 + v1(i)**2) / RG & + 0.11 * 14e-6 / SQRT(cdragm(i) * tmp) z0m(i) = MAX(1.5e-05, z0m(i)) z0h(i) = 0.4 * 14e-6 / SQRT(cdragm(i) * tmp) ENDDO ELSE IF (iflag_z0_oce==-1) THEN DO i = 1, knon z0m(i) = z0min z0h(i) = z0min ENDDO ELSE CALL abort_physic(modname, 'version non prevue', 1) ENDIF IF (activate_ocean_skin >= 1) THEN IF (type_ocean /= 'couple') sss(:knon) = 35. CALL bulk_flux(tkt, tks, taur, dter, dser, t_int, s_int, ds_ns, dt_ns, & u = windsp(:knon), t_ocean_1 = tsurf_new(:knon), s1 = sss(:knon), & rain = precip_rain(:knon) + precip_totsnow(:knon), & hf = - fluxsens(:knon), hlb = - fluxlat(:knon), & rnl = - lwnet(:knon), & tau = sqrt(flux_u1(:knon)**2 + flux_v1(:knon)**2), rhoa = rhoa, & xlv = [(rlvtt, i = 1, knon)], rf = - sens_prec_liq, dtime = dtime, & rns = swnet(:knon)) delta_sst = t_int - tsurf_new(:knon) delta_sal = s_int - sss(:knon) IF (activate_ocean_skin == 2) THEN tsurf_new(:knon) = t_int IF (type_ocean == 'couple') dt_ds = (tks / tkt) * dter end if end if END SUBROUTINE surf_ocean !**************************************************************************** END MODULE surf_ocean_mod