! ! $Id: surf_ocean_mod.F90 5662 2025-05-20 14:24:41Z fairhead $ ! 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, & !GG dt_ds, tkt, tks, taur, sss) dt_ds, tkt, tks, taur, sss, & dthetadz300,Ampl & !GG #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 #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 clesphys_mod_h USE yomcst_mod_h USE limit_read_mod USE config_ocean_skin_m, ONLY: activate_ocean_skin ! ! 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. ! for cycle_diurne and for iflag_z0_oce==-1 (prescribed z0) ! 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 !GG REAL, DIMENSION(klon), INTENT(IN) :: dthetadz300 REAL, DIMENSION(klon), INTENT(OUT) :: Ampl ! ! 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 .ne. 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, & !GG tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa) tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa, & dthetadz300, pctsrf, Ampl & !GG #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.NE.'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.GE.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