! ! $Id: ocean_forced_mod.F90 5662 2025-05-20 14:24:41Z fairhead $ ! MODULE ocean_forced_mod ! ! This module is used for both the sub-surfaces ocean and sea-ice for the case of a ! forced ocean, "ocean=force". ! IMPLICIT NONE CONTAINS ! !**************************************************************************************** ! SUBROUTINE ocean_forced_noice( & itime, dtime, jour, knon, knindex, & p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, & 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, & dthetadz300,pctsrf,Ampl & #ifdef ISO ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, & xtsnow,xtevap,h1 & #endif ) ! ! This subroutine treats the "open ocean", all grid points that are not entierly covered ! by ice. ! The routine receives data from climatologie file limit.nc and does some calculations at the ! surface. ! USE dimphy USE calcul_fluxs_mod USE limit_read_mod USE mod_grid_phy_lmdz USE indice_sol_mod USE surface_data, ONLY : iflag_leads USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o use config_ocean_skin_m, only: activate_ocean_skin #ifdef ISO USE infotrac_phy, ONLY: ntiso,niso USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall #ifdef ISOVERIF USE isotopes_mod, ONLY: iso_eau,ridicule !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix USE isotopes_verif_mod #endif #endif USE flux_arp_mod_h USE clesphys_mod_h USE yomcst_mod_h ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: itime, jour, knon INTEGER, DIMENSION(klon), INTENT(IN) :: knindex REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon), INTENT(IN) :: p1lay REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow 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) :: tsurf_in real, intent(in):: rhoa(:) ! (knon) density of moist air (kg / m3) !GG REAL, DIMENSION(klon), INTENT(IN) :: dthetadz300 REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf ! #ifdef ISO REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum REAL, DIMENSION(klon), INTENT(IN) :: rlat #endif ! In/Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(INOUT) :: radsol REAL, DIMENSION(klon), INTENT(INOUT) :: snow REAL, DIMENSION(klon), INTENT(INOUT) :: agesno !? put to 0 in ocean #ifdef ISO REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsnow REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce #endif ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: qsurf REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l REAL, INTENT(out):: sens_prec_liq(:) ! (knon) !GG REAL, DIMENSION(klon), INTENT(OUT) :: Ampl ! #ifdef ISO REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux REAL, DIMENSION(klon), INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation #endif ! Local variables !**************************************************************************************** INTEGER :: i, j REAL, DIMENSION(klon) :: cal, beta, dif_grnd REAL, DIMENSION(klon) :: alb_neig, tsurf_lim, zx_sl REAL, DIMENSION(klon) :: u0, v0 REAL, DIMENSION(klon) :: u1_lay, v1_lay LOGICAL :: check=.FALSE. REAL, DIMENSION(knon) :: sens_prec_sol REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol ! GG REAL, DIMENSION(klon) :: l_CBL, sicfra ! #ifdef ISO REAL, PARAMETER :: t_coup = 273.15 #endif !**************************************************************************************** ! Start calculation !**************************************************************************************** IF (check) WRITE(*,*)' Entering ocean_forced_noice' #ifdef ISO #ifdef ISOVERIF DO i = 1, knon IF (iso_eau > 0) THEN CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), & & spechum(i),'ocean_forced_mod 111', & & errmax,errmaxrel) CALL iso_verif_egalite_choix(snow(i), & & xtsnow(iso_eau,i),'ocean_forced_mod 117', & & errmax,errmaxrel) ENDIF !IF (iso_eau > 0) THEN ENDDO !DO i=1,knon #endif #endif !**************************************************************************************** ! 1) ! Read sea-surface temperature from file limit.nc ! !**************************************************************************************** !--sb: !!jyg if (knon.eq.1) then ! single-column model if (klon_glo.eq.1) then ! single-column model ! EV: now surface Tin flux_arp.h !CALL read_tsurf1d(knon,tsurf_lim) ! new DO i = 1, knon tsurf_lim(i) = tg ENDDO else ! GCM CALL limit_read_sst(knon,knindex,tsurf_lim & #ifdef ISO & ,Roce,rlat & #endif & ) endif ! knon !sb-- !**************************************************************************************** ! 2) ! Flux calculation ! !**************************************************************************************** ! Set some variables for calcul_fluxs !cal = 0. !beta = 1. !dif_grnd = 0. ! EV: use calbeta to calculate beta ! Need to initialize qsurf for calbeta but it is not modified by this routine qsurf(:)=0. CALL calbeta(dtime, is_oce, knon, snow, qsurf, beta, cal, dif_grnd) alb_neig(:) = 0. agesno(:) = 0. lat_prec_liq = 0.; lat_prec_sol = 0. ! Suppose zero surface speed u0(:)=0.0 v0(:)=0.0 u1_lay(:) = u1(:) - u0(:) v1_lay(:) = v1(:) - v0(:) ! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf CALL calcul_fluxs(knon, is_oce, dtime, & merge(tsurf_in, tsurf_lim, activate_ocean_skin == 2), p1lay, cal, & beta, cdragh, cdragq, ps, & precip_rain, precip_snow, snow, qsurf, & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) if (activate_ocean_skin == 2) tsurf_new = tsurf_lim do j = 1, knon i = knindex(j) sens_prec_liq_o(i,1) = sens_prec_liq(j) sens_prec_sol_o(i,1) = sens_prec_sol(j) lat_prec_liq_o(i,1) = lat_prec_liq(j) lat_prec_sol_o(i,1) = lat_prec_sol(j) enddo !GG if (iflag_leads == 1) then l_CBL = -52381.*dthetadz300 + 3008.1 Ampl = 6.012e-08*l_CBL**2 - 4.036e-04*l_CBL + 1.4979 WHERE(Ampl(:)>1.2) Ampl(:)=1.2 sicfra(:)=pctsrf(:,is_sic)/(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter)) WHERE(pctsrf(:,is_sic)+pctsrf(:,is_oce)<EPSFRA) sicfra(:)=0. WHERE(sicfra<0.7) Ampl(:)=1. WHERE((sicfra>0.7).and.(sicfra<0.9)) Ampl=((sicfra-0.7)/0.2)*Ampl+((0.9-sicfra)/0.2) fluxsens=Ampl*fluxsens dflux_s=Ampl*dflux_s endif ! - Flux calculation at first modele level for U and V CALL calcul_flux_wind(knon, dtime, & u0, v0, u1, v1, gustiness, cdragm, & AcoefU, AcoefV, BcoefU, BcoefV, & p1lay, temp_air, & flux_u1, flux_v1) #ifdef ISO CALL calcul_iso_surf_oce_vectall(klon, knon,t_coup, & & ps,tsurf_new,spechum,u1_lay, v1_lay, xtspechum, & & evap, Roce,xtevap,h1 & #ifdef ISOTRAC & ,knindex & #endif & ) #endif #ifdef ISO #ifdef ISOVERIF ! write(*,*) 'ocean_forced_mod 176: sortie de ocean_forced_noice' IF (iso_eau > 0) THEN DO i = 1, knon CALL iso_verif_egalite_choix(snow(i), & & xtsnow(iso_eau,i),'ocean_forced_mod 180', & & errmax,errmaxrel) ENDDO ! DO j=1,knon ENDIF !IF (iso_eau > 0) THEN #endif #endif END SUBROUTINE ocean_forced_noice ! !*************************************************************************************** ! SUBROUTINE ocean_forced_ice( & itime, dtime, jour, knon, knindex, & tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air,spechum, & AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & !GG ps, u1, v1, gustiness, & ps, u1, v1, gustiness, pctsrf, & !GG radsol, snow, qsol, agesno, tsoil, & qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & !GG tsurf_new, dflux_s, dflux_l, rhoa) tsurf_new, dflux_s, dflux_l, rhoa, swnet, hice, tice, bilg_cumul, & fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, & dtice_melt, dtice_snow2sic & !GG #ifdef ISO ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, & xtsnow, xtsol,xtevap,Rland_ice & #endif ) ! ! This subroutine treats the ocean where there is ice. ! The routine reads data from climatologie file and does flux calculations at the ! surface. ! USE dimphy USE geometry_mod, ONLY: longitude,latitude USE calcul_fluxs_mod !GG USE surface_data, ONLY : calice, calsno USE surface_data, ONLY : calice, calsno, iflag_seaice, iflag_seaice_alb, & sice_cond, sisno_cond, sisno_den, sisno_min, sithick_min, sisno_wfact, & amax_s,amax_n, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, & si_pen_frac, si_pen_ext, fseaN, fseaS, iflag_leads USE geometry_mod, ONLY: longitude,latitude,latitude_deg !GG USE limit_read_mod USE fonte_neige_mod, ONLY : fonte_neige USE indice_sol_mod USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o #ifdef ISO USE infotrac_phy, ONLY: niso, ntiso USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall #ifdef ISOVERIF USE isotopes_mod, ONLY: iso_eau,ridicule !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix USE isotopes_verif_mod #endif #endif USE flux_arp_mod_h USE clesphys_mod_h USE yomcst_mod_h USE dimsoil_mod_h, ONLY: nsoilmx ! INCLUDE "indicesol.h" ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: itime, jour, knon INTEGER, DIMENSION(klon), INTENT(IN) :: knindex REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in REAL, DIMENSION(klon), INTENT(IN) :: p1lay REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow 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, intent(in):: rhoa(:) ! (knon) density of moist air (kg / m3) !GG REAL, DIMENSION(klon), INTENT(IN) :: swnet REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf !GG #ifdef ISO REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum REAL, DIMENSION(niso,klon), INTENT(IN) :: Roce REAL, DIMENSION(niso,klon), INTENT(IN) :: Rland_ice #endif ! In/Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(INOUT) :: radsol REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol REAL, DIMENSION(klon), INTENT(INOUT) :: agesno REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil !GG REAL, DIMENSION(klon), INTENT(INOUT) :: hice REAL, DIMENSION(klon), INTENT(INOUT) :: tice REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul REAL, DIMENSION(klon), INTENT(INOUT) :: fcds REAL, DIMENSION(klon), INTENT(INOUT) :: fcdi REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_growth REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_melt REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt REAL, DIMENSION(klon), INTENT(INOUT) :: dh_snow2sic REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_melt REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_snow2sic !GG #ifdef ISO REAL, DIMENSION(niso,klon), INTENT(INOUT) :: xtsnow REAL, DIMENSION(niso,klon), INTENT(IN) :: xtsol #endif ! Output arguments !**************************************************************************************** REAL, DIMENSION(klon), INTENT(OUT) :: qsurf 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(klon), INTENT(OUT) :: evap, fluxsens, fluxlat REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l #ifdef ISO REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap #endif ! Local variables !**************************************************************************************** LOGICAL :: check=.FALSE. INTEGER :: i, j REAL :: zfra REAL, PARAMETER :: t_grnd=271.35 REAL, DIMENSION(klon) :: cal, beta, dif_grnd, capsol, icesub REAL, DIMENSION(klon) :: alb_neig, tsurf_tmp REAL, DIMENSION(klon) :: soilcap, soilflux REAL, DIMENSION(klon) :: u0, v0 REAL, DIMENSION(klon) :: u1_lay, v1_lay REAL, DIMENSION(knon) :: sens_prec_liq, sens_prec_sol REAL, DIMENSION(klon) :: lat_prec_liq, lat_prec_sol !GG INTEGER :: ki INTEGER :: cpl_pas REAL, DIMENSION(klon) :: bilg, fsic, f_bot REAL, PARAMETER :: latent_ice = 334.0e3 REAL, PARAMETER :: rau_ice = 917.0 REAL, PARAMETER :: kice=2.2 REAL :: f_cond, f_swpen, f_cond_s, f_cond_i REAL :: ustar, uscap, ustau ! for snow/ice albedo: REAL :: alb_snow, alb_ice, alb_pond REAL :: frac_snow, frac_ice, frac_pond REAL :: z1_i, z2_i, z1_s, zlog ! height parameters ! for ice melt / freeze REAL :: e_melt, snow_evap, h_test ! dhsic, dfsic change in ice mass, fraction. REAL :: dhsic, dfsic, frac_mf REAL :: fsea, amax REAL :: hice_i, tice_i, fsic_new ! snow and ice physical characteristics: REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp REAL, PARAMETER :: t_melt=273.15 ! melting ice temp REAL :: sno_den!=sisno_den !mean snow density, kg/m3 REAL, PARAMETER :: ice_den=917. ! ice density REAL, PARAMETER :: sea_den=1025. ! sea water density REAL :: ice_cond!=sice_cond*ice_den !conductivity of ice REAL :: sno_cond!=sisno_cond*sno_den ! conductivity of snow REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice REAL, PARAMETER :: sea_cap=3995. ! specific heat capacity, water REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice ! control of snow and ice cover & freeze / melt (heights converted to kg/m2) REAL :: snow_min!=sisno_min*sno_den !critical snow height 5 cm REAL :: snow_wfact!=sisno_wfact ! max fraction of falling snow blown into ocean REAL, PARAMETER :: ice_frac_min=0.005 REAL :: h_ice_min!=sithick_min ! min ice thickness ! below ice_thin, priority is melt lateral / grow height ! ice_thin is also height of new ice REAL, PARAMETER :: h_ice_max=7 ! max ice height ! Ice thickness parameter for lateral growth REAL, PARAMETER :: h_ice_thick=1.5 REAL, PARAMETER :: h_ice_thin=0.15 ! albedo and radiation parameters INTEGER, SAVE :: iflag_sic_albedo ! albedo old or NEMO REAL :: alb_sno_dry!=rn_alb_sdry !dry snow albedo REAL :: alb_sno_wet!=rn_alb_smlt !wet snow albedo REAL :: alb_ice_dry!=rn_alb_idry !dry thick ice REAL :: alb_ice_wet!=rn_alb_imlt !melting thick ice ! new (Toyoda 2020) albedo ! Values for snow / ice, dry / melting, visible / near IR REAL, PARAMETER :: alb_sdry_vis=0.98 REAL, PARAMETER :: alb_smlt_vis=0.88 REAL, PARAMETER :: alb_sdry_nir=0.7 REAL, PARAMETER :: alb_smlt_nir=0.55 REAL, PARAMETER :: alb_idry_vis=0.78 REAL, PARAMETER :: alb_imlt_vis=0.705 REAL, PARAMETER :: alb_idry_nir=0.36 REAL, PARAMETER :: alb_imlt_nir=0.285 REAL, PARAMETER :: h_ice_alb=0.5*ice_den ! height for full ice albedo REAL, PARAMETER :: h_sno_alb=0.02*300 ! height for control of snow fraction REAL :: pen_frac !=si_pen_frac !fraction of shortwave penetrating into the ! ice (not snow). Should be visible only, not NIR REAL :: pen_ext !=si_pen_ext !extinction length of penetrating shortwave (m-1) ! HF from ocean below ice ! REAL, PARAMETER :: fseaN=2.0 ! NH ! REAL, PARAMETER :: fseaS=4.0 ! SH !GG #ifdef ISO REAL, PARAMETER :: t_coup = 273.15 REAL, DIMENSION(klon) :: fq_fonte_diag REAL, DIMENSION(klon) :: fqfonte_diag REAL, DIMENSION(klon) :: snow_evap_diag REAL, DIMENSION(klon) :: fqcalving_diag REAL, DIMENSION(klon) :: run_off_lic_diag REAL :: coeff_rel_diag REAL :: max_eau_sol_diag REAL, DIMENSION(klon) :: runoff_diag INTEGER IXT REAL, DIMENSION(niso,klon) :: xtsnow_prec, xtsol_prec REAL, DIMENSION(klon) :: snow_prec, qsol_prec #endif !**************************************************************************************** ! Start calculation !**************************************************************************************** IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon !**************************************************************************************** ! 1) ! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, flux_u1, flux_v1 ! dflux_s, dflux_l and qsurf !**************************************************************************************** tsurf_tmp(:) = tsurf_in(:) !GG IF (iflag_seaice==0) THEN ! Old LMDZ sea ice surface !GG ! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd) IF (soil_model) THEN ! update tsoil and calculate soilcap and soilflux CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, & & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil,soilcap, soilflux) cal(1:knon) = RCPD / soilcap(1:knon) radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) dif_grnd = 1.0 / tau_gl ELSE dif_grnd = 1.0 / tau_gl cal = RCPD * calice WHERE (snow > 0.0) cal = RCPD * calsno ENDIF !GG ELSEIF (iflag_seaice==2) THEN sno_den=sisno_den !mean snow density, kg/m3 ice_cond=sice_cond*ice_den !conductivity of ice sno_cond=sisno_cond*sno_den ! conductivity of snow snow_min=sisno_min*sno_den !critical snow height 5 cm snow_wfact=sisno_wfact ! max fraction of falling snow blown into ocean h_ice_min=sithick_min ! min ice thickness alb_sno_dry=rn_alb_sdry !dry snow albedo alb_sno_wet=rn_alb_smlt !wet snow albedo alb_ice_dry=rn_alb_idry !dry thick ice alb_ice_wet=rn_alb_imlt !melting thick ice pen_frac=si_pen_frac !fraction of shortwave penetrating into the pen_ext=si_pen_ext !extinction length of penetrating shortwave (m-1) bilg(:)=0. dif_grnd(:)=0. beta(:) = 1. fsic(:) = pctsrf(:,is_sic) cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour ! Surface, snow-ice and ice-ocean fluxes. ! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd) !write(*,*) 'radsol 1',radsol(1:100) DO i=1,knon ki=knindex(i) IF (snow(i).GT.snow_min) THEN ! 1 / snow-layer heat capacity cal(i)=2.*RCPD/(snow(i)*ice_cap) ! adjustment time-scale of conductive flux dif_grnd(i) = cal(i) * sno_cond / snow(i) / RCPD ! for conductive flux f_cond_s = sno_cond * (tice(ki)-t_freeze) / snow(i) radsol(i) = radsol(i)+f_cond_s ! all shortwave flux absorbed f_swpen=0. ELSE ! bare ice. f_cond_s = 0. ! 1 / ice-layer heat capacity cal(i) = 2.*RCPD/(hice(ki)*ice_den*ice_cap) ! adjustment time-scale of conductive flux dif_grnd(i) = cal(i) * ice_cond / (hice(ki)*ice_den) / RCPD ! penetrative shortwave flux... f_swpen=swnet(i)*pen_frac*exp(-pen_ext*hice(ki)) radsol(i) = radsol(i)-f_swpen ! GG no conductive flux in this case? END IF bilg(ki)=f_swpen-f_cond_s END DO endif ! beta = 1.0 lat_prec_liq = 0.; lat_prec_sol = 0. ! Suppose zero surface speed u0(:)=0.0 v0(:)=0.0 u1_lay(:) = u1(:) - u0(:) v1_lay(:) = v1(:) - v0(:) CALL calcul_fluxs(knon, is_sic, dtime, & tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, & precip_rain, precip_snow, snow, qsurf, & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) do j = 1, knon i = knindex(j) sens_prec_liq_o(i,2) = sens_prec_liq(j) sens_prec_sol_o(i,2) = sens_prec_sol(j) lat_prec_liq_o(i,2) = lat_prec_liq(j) lat_prec_sol_o(i,2) = lat_prec_sol(j) enddo ! - Flux calculation at first modele level for U and V CALL calcul_flux_wind(knon, dtime, & u0, v0, u1, v1, gustiness, cdragm, & AcoefU, AcoefV, BcoefU, BcoefV, & p1lay, temp_air, & flux_u1, flux_v1) !**************************************************************************************** ! 2) ! Calculations due to snow and runoff ! !**************************************************************************************** !GG if (iflag_seaice==0) then !GG #ifdef ISO ! verif #ifdef ISOVERIF DO i = 1, knon IF (iso_eau > 0) THEN IF (snow(i) > ridicule) THEN CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), & & 'interfsurf 964',errmax,errmaxrel) ENDIF !IF ((snow(i) > ridicule)) THEN ENDIF !IF (iso_eau > 0) THEN ENDDO !DO i=1,knon #endif ! end verif DO i = 1, knon snow_prec(i) = snow(i) DO ixt = 1, niso xtsnow_prec(ixt,i) = xtsnow(ixt,i) ENDDO !DO ixt=1,niso ! initialisation: fq_fonte_diag(i) = 0.0 fqfonte_diag(i) = 0.0 snow_evap_diag(i)= 0.0 ENDDO !DO i=1,knon #endif CALL fonte_neige( knon, is_sic, knindex, dtime, & tsurf_tmp, precip_rain, precip_snow, & snow, qsol, tsurf_new, evap, icesub & #ifdef ISO & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag & & ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag & #endif & ) #ifdef ISO ! isotopes: tout est externalisé !#ifdef ISOVERIF ! write(*,*) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall' ! write(*,*) 'klon,knon=',klon,knon !#endif CALL calcul_iso_surf_sic_vectall(klon,knon, & & evap,snow_evap_diag,Tsurf_new,Roce,snow, & & fq_fonte_diag,fqfonte_diag,dtime,t_coup, & & precip_snow,xtprecip_snow,xtprecip_rain, snow_prec,xtsnow_prec, & & xtspechum,spechum,ps, & & xtevap,xtsnow,fqcalving_diag, & & knindex,is_sic,run_off_lic_diag,coeff_rel_diag,Rland_ice & & ) #ifdef ISOVERIF !write(*,*) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall' IF (iso_eau > 0) THEN DO i = 1, knon CALL iso_verif_egalite_choix(snow(i), & & xtsnow(iso_eau,i),'ocean_forced_mod 396', & & errmax,errmaxrel) ENDDO ! DO j=1,knon ENDIF !IF (iso_eau > 0) then #endif !#ifdef ISOVERIF #endif !#ifdef ISO ! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno) ! CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0. alb1_new(:) = 0.0 DO i=1, knon zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0))) alb1_new(i) = alb_neig(i) * zfra + 0.6 * (1.0-zfra) ENDDO alb2_new(:) = alb1_new(:) !GG else DO i=1,knon ki=knindex(i) ! ocean-ice heat flux fsea=fseaS amax=amax_s if (latitude(ki)>0) THEN fsea=fseaN amax=amax_n ENDIF IF (snow(i).GT.snow_min) THEN ! snow conductive flux after calcul_fluxs (pos up) f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i) ! 1 / heat capacity and conductive timescale uscap = 2. / ice_cap / (snow(i)+hice(ki)*ice_den) ustau = uscap * ice_cond / (hice(ki)*ice_den) ! update ice temp tice(ki) = (tice(ki) + dtime*(ustau*t_freeze - uscap*f_cond_s)) & / (1 + dtime*ustau) ELSE ! bare ice tice(ki)=tsurf_new(i) ENDIF ! ice conductive flux (pos up) f_cond_i = ice_cond * (t_freeze-tice(ki)) / (hice(ki)*ice_den) f_bot(i) = fsea - f_cond_i fcdi(ki) = f_cond_i - fsea fcds(i) = f_cond_s !bilg(ki) = bilg(ki)+f_cond_i END DO !**************************************************************************************** ! 2) Update snow and ice surface : thickness !**************************************************************************************** IF (iflag_seaice==1) THEN ! Read from limit CALL limit_read_hice(knon,knindex,hice) ENDIF ! Formula Krinner et al. 1997 : h = (0.2 + 3.8(f_min**2))(1 + 2(f- f_min)) DO i=1,knon ki=knindex(i) IF (precip_snow(i) > 0.) THEN snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki))) END IF ! snow and ice sublimation IF (evap(i) > 0.) THEN snow_evap = MIN (snow(i) / dtime, evap(i)) snow(i) = snow(i) - snow_evap * dtime snow(i) = MAX(0.0, snow(i)) IF (iflag_seaice==2) THEN hice(ki) = MAX(0.0,hice(ki)-(evap(i)-snow_evap)*dtime/ice_den) ENDIF ENDIF ! Melt / Freeze snow from above if Tsurf>0 IF (tsurf_new(i).GT.t_melt) THEN ! energy available for melting snow (in kg of melted snow /m2) e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. & /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i)) ! remove snow tice_i=tice(ki) IF (snow(i).GT.e_melt) THEN snow(i)=snow(i)-e_melt tsurf_new(i)=t_melt ELSE ! all snow is melted ! add remaining heat flux to ice e_melt=e_melt-snow(i) tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*hice(ki)*ice_den) tsurf_new(i)=tice(ki) END IF dtice_melt(ki)=(tice(ki)-tice_i)/dtime END IF ! Bottom melt / grow ! bottom freeze if bottom flux (cond + oce-ice) <0 IF (iflag_seaice==2) THEN IF (f_bot(i).LT.0) THEN ! larger fraction of bottom growth frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thick) & / (h_ice_max-h_ice_thick))) ! quantity of new ice (formed at mean ice temp) e_melt= -f_bot(i) * dtime * fsic(ki) & / (ice_lat+ice_cap/2.*(t_freeze-tice(ki))) ! first increase height to h_thick dhsic=MAX(0.,MIN(h_ice_thick-hice(ki),e_melt/(fsic(ki)*ice_den))) hice_i=hice(ki) hice(ki)=dhsic+hice(ki) e_melt=e_melt-fsic(ki)*dhsic IF (e_melt.GT.0.) THEN ! frac_mf fraction used for lateral increase dfsic=MIN(amax-fsic(ki),e_melt*frac_mf/ (hice(ki)*ice_den) ) ! No lateral growth -> forced ocean !fsic(ki)=fsic(ki)+dfsic e_melt=e_melt-dfsic*hice(ki)*ice_den ! rest used to increase height hice(ki)=MIN(h_ice_max,hice(ki)+e_melt/( fsic(ki) * ice_den ) ) END IF dh_basal_growth(ki)=(hice(ki)-hice_i)/dtime ! melt from below if bottom flux >0 ELSE ! larger fraction of lateral melt from warm ocean frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thin) & / (h_ice_thick-h_ice_thin))) ! bring ice to freezing and melt from below ! quantity of melted ice e_melt= f_bot(i) * dtime * fsic(ki) & / (ice_lat+ice_cap/2.*(tice(ki)-t_freeze)) ! first decrease height to h_thick hice_i=hice(ki) dhsic=MAX(0.,MIN(hice(ki)-h_ice_thick,e_melt/(fsic(ki)*ice_den))) hice(ki)=hice(ki)-dhsic e_melt=e_melt-fsic(ki)*dhsic*ice_den IF (e_melt.GT.0) THEN ! frac_mf fraction used for height decrease dhsic=MAX(0.,MIN(hice(ki)-h_ice_min,e_melt/ice_den*frac_mf/fsic(ki))) hice(ki)=hice(ki)-dhsic e_melt=e_melt-fsic(ki)*dhsic*ice_den ! rest used to decrease fraction (up to 0!) dfsic=MIN(fsic(ki),e_melt/(hice(ki)*ice_den)) ! Remaining heat not used if everything melted e_melt=e_melt-dfsic*hice(ki)*ice_den ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime END IF dh_basal_melt(ki)=(hice(ki)-hice_i)/dtime END IF END IF ! melt ice from above if Tice>0 tice_i=tice(ki) IF (tice(ki).GT.t_melt) THEN IF (iflag_seaice==2) THEN ! quantity of ice melted (kg/m2) e_melt=MAX(hice(ki)*ice_den*(tice(ki)-t_melt)*ice_cap/2. & /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0) ! melt from above, height only hice_i=hice(ki) dhsic=MIN(hice(ki)-h_ice_min,e_melt/ice_den) dh_top_melt(i)=dhsic e_melt=e_melt-dhsic IF (e_melt.GT.0) THEN ! lateral melt if ice too thin dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/(h_ice_min*ice_den)*fsic(ki)) ! if all melted do nothing with remaining heat e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min*ice_den) ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime END IF hice(ki)=hice(ki)-dhsic dh_top_melt(ki)=(hice(ki)-hice_i)/dtime ! surface temperature at melting point END IF tice(ki)=t_melt tsurf_new(i)=t_melt END IF dtice_melt(ki)=dtice_melt(ki)+tice(ki)-tice_i ! convert snow to ice if below floating line h_test=(hice(ki)*ice_den+snow(i))-hice(ki)*sea_den IF ((h_test.GT.0.).AND.(hice(ki).GT.h_ice_min)) THEN !snow under water ! extra snow converted to ice (with added frozen sea water) IF (iflag_seaice==2) THEN dhsic=h_test/(sea_den-ice_den+sno_den) hice(ki)=hice(ki)+dhsic ENDIF snow(i)=snow(i)-dhsic*sno_den ! available energy (freeze sea water + bring to tice) e_melt=dhsic*ice_den*(1.-sno_den/ice_den)*(ice_lat+ & ice_cap/2.*(t_freeze-tice(ki))) ! update ice temperature tice_i=tice(ki) tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+hice(ki)*ice_den) IF (iflag_seaice==2) THEN dh_snow2sic(ki)=dhsic/dtime END IF dtice_snow2sic(ki)=(tice(ki)-tice_i)/dtime END IF END DO !write(*,*) 'hice 2',hice(1:100) !write(*,*) 'tice 2',tice(1:100) iflag_sic_albedo=iflag_seaice_alb !******************************************************************************* ! 3) cumulate ice-ocean fluxes, update tslab, lateral grow !***********************************************o******************************* !cumul fluxes. bilg_cumul(:)=bilg_cumul(:)+bilg(:)/float(cpl_pas) IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab bilg_cumul(:)=0. END IF ! coupling time ! write(*,*) 'hice 3',hice(1:100) ! write(*,*) 'tice 3',tice(1:100) !tests ice fraction WHERE (fsic.LT.ice_frac_min) tice=t_melt hice=h_ice_min END WHERE !write(*,*) 'hice 4',hice(1:100) !write(*,*) 'tice 4',tice(1:100) endif !**************************************************************************************** ! 4) Compute sea-ice and snow albedo !**************************************************************************************** IF (iflag_seaice > 0) THEN SELECT CASE (iflag_sic_albedo) CASE(0) ! old slab albedo : single value. age of snow, melt ponds. DO i=1,knon ki=knindex(i) ! snow albedo: update snow age IF (snow(i).GT.0.0001) THEN agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)& * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.) ELSE agesno(i)=0.0 END IF ! snow albedo alb_snow=alb_sno_wet+(alb_sno_dry-alb_sno_wet)*EXP(-agesno(i)/50.) ! ice albedo (varies with ice tkickness and temp) alb_ice=MAX(0.0,0.13*LOG(100.*hice(ki))+0.1) !alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1) IF (tice(ki).GT.t_freeze-0.01) THEN alb_ice=MIN(alb_ice,alb_ice_wet) ELSE alb_ice=MIN(alb_ice,alb_ice_dry) END IF ! pond albedo alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) ! pond fraction frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) ! snow fraction frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min)) ! ice fraction frac_ice=MAX(0.0,1.-frac_pond-frac_snow) ! total albedo alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond END DO alb2_new(:) = alb1_new(:) CASE(1) ! New visible and IR albedos, dry / melting snow ! based on Toyoda et al, 2020 DO i=1,knon ki=knindex(i) ! snow fraction frac_snow = snow(i) / (snow(i) + h_sno_alb) ! dependence of ice albedo with ice thickness frac_ice = MIN(1.,ATAN(4.*hice(ki)*ice_den) / ATAN(4.*h_ice_alb)) ! Total (for ice, min = 0.066 = alb_ocean) IF (tice(ki).GT.t_melt) THEN alb_ice = 0.066 + (alb_imlt_vis - 0.066)*frac_ice alb1_new(i)=alb_smlt_vis*frac_snow + alb_ice*(1.-frac_snow) alb_ice = 0.066 + (alb_imlt_nir - 0.066)*frac_ice alb2_new(i)=alb_smlt_nir*frac_snow + alb_ice*(1.-frac_snow) ELSEIF (tice(ki).GT.t_melt - 1.) THEN frac_pond = tice(ki) - t_freeze alb_snow = alb_smlt_vis*frac_pond + alb_sdry_vis*(1.-frac_pond) alb_ice = alb_imlt_vis*frac_pond + alb_idry_vis*(1.-frac_pond) alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice alb1_new(i)= alb_snow*frac_snow + alb_ice*(1.-frac_snow) alb_snow = alb_smlt_nir*frac_pond + alb_sdry_nir*(1.-frac_pond) alb_ice = alb_imlt_nir*frac_pond + alb_idry_nir*(1.-frac_pond) alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice alb2_new(i)= alb_snow*frac_snow + alb_ice*(1.-frac_snow) ELSE alb_ice = 0.066 + (alb_idry_vis - 0.066)*frac_ice alb1_new(i)=alb_sdry_vis*frac_snow + alb_ice*(1.-frac_snow) alb_ice = 0.066 + (alb_idry_nir - 0.066)*frac_ice alb2_new(i)=alb_sdry_nir*frac_snow + alb_ice*(1.-frac_snow) ENDIF END DO CASE(2) ! LIM3 scheme. Uses clear sky / overcast value, with 50% clear sky z1_i = 1.5 * ice_den z2_i = 0.05 * ice_den zlog = 1. / (LOG(1.5) - LOG(0.05)) z1_s = 1. / (0.025 * sno_den) DO i=1,knon ki=knindex(i) ! temperature above / below 0 IF (tice(ki).GE.t_melt) THEN alb_ice = alb_ice_wet alb_snow = alb_sno_wet ELSE alb_ice = alb_ice_dry alb_snow = alb_sno_dry ENDIF ! ice thickness IF (hice(ki)*ice_den.LT.z2_i) THEN alb_ice = 0.066 + 0.114 * hice(ki)*ice_den / z2_i ELSEIF (hice(ki)*ice_den.LT.z1_i) THEN alb_ice = alb_ice + (0.18 - alb_ice) & * (LOG(z1_i) - LOG(hice(ki)*ice_den)) * zlog ENDIF ! ice or snow depending on snow thickness alb_snow = alb_snow - (alb_snow -alb_ice) * EXP(- snow(i) * z1_s) ! Effect of clouds (polynomial fit with 50% clouds) alb1_new(i) = alb_snow - 0.5 * (-0.1010 * alb_snow*alb_snow & + 0.1933*alb_snow - 0.0148) alb2_new(i) = alb1_new(i) END DO CASE(3) CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0. alb1_new(:) = 0.0 DO i=1, knon zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0))) alb1_new(i) = alb_neig(i) * zfra + 0.6 * (1.0-zfra) ENDDO alb2_new(:) = alb1_new(:) print*,'alb_neig=',alb_neig print*,'zfra=',zfra print*,'snow=',snow print*,'alb1_new=',alb1_new print*,'alb2_new=',alb2_new END SELECT END IF ! ------ End Albedo ---------- !GG END SUBROUTINE ocean_forced_ice !************************************************************************ ! 1D case !************************************************************************ ! SUBROUTINE read_tsurf1d(knon,sst_out) ! ! This subroutine specifies the surface temperature to be used in 1D simulations ! ! USE dimphy, ONLY : klon ! ! INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid ! REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model ! ! INTEGER :: i ! COMMON defined in lmdz1d.F: ! real ts_cur ! common /sst_forcing/ts_cur ! ! DO i = 1, knon ! sst_out(i) = ts_cur ! ENDDO ! ! END SUBROUTINE read_tsurf1d ! ! !************************************************************************ END MODULE ocean_forced_mod