- Timestamp:
- Feb 6, 2026, 10:58:17 AM (2 months ago)
- Location:
- LMDZ6/trunk/libf/phylmd
- Files:
-
- 5 edited
-
conf_phys_m.f90 (modified) (6 diffs)
-
ocean_forced_mod.F90 (modified) (7 diffs)
-
phyetat0_mod.f90 (modified) (2 diffs)
-
surf_ocean_mod.F90 (modified) (2 diffs)
-
surface_data.f90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/conf_phys_m.f90
r6059 r6070 199 199 INTEGER, SAVE :: iflag_cycle_diurne_omp 200 200 LOGICAL, SAVE :: soil_model_omp,liqice_in_radocond_omp 201 !GG202 201 INTEGER,SAVE :: iflag_seaice_omp, iflag_seaice_alb_omp 203 202 INTEGER,SAVE :: iflag_leads_omp … … 208 207 REAL,SAVE :: si_pen_frac_omp,si_pen_ext_omp 209 208 REAL,SAVE :: fseaN_omp, fseaS_omp 210 !GG209 REAL,SAVE :: sic_hice_fixed_omp 211 210 LOGICAL, SAVE :: ok_limitvrai_omp 212 211 INTEGER, SAVE :: nbapp_rad_omp, iflag_con_omp … … 902 901 CALL getin('iflag_seaice_alb',iflag_seaice_alb_omp) 903 902 903 !Config Key = sic_hice_fixed 904 !Config Desc = imposed sea-ice thickness [m] 905 !Config Def = 906 sic_hice_fixed_omp = 1.0 907 CALL getin('sic_hice_fixed',sic_hice_fixed_omp) 908 909 904 910 !Config Key = iflag_seaice_alb 905 911 !Config Desc = Flag de sea ice leads … … 2179 2185 soil_model = soil_model_omp 2180 2186 liqice_in_radocond = liqice_in_radocond_omp 2181 ! GG 2187 sic_hice_fixed = sic_hice_fixed_omp 2182 2188 sice_cond = sice_cond_omp 2183 2189 sisno_cond = sisno_cond_omp … … 2710 2716 WRITE(lunout,*) ' var_fco2_land_cor = ', var_fco2_land_cor 2711 2717 WRITE(lunout,*) ' dms_cycle_cpl = ', dms_cycle_cpl 2712 !GG2713 2718 WRITE(lunout,*) ' iflag_seaice = ', iflag_seaice 2714 2719 WRITE(lunout,*) ' iflag_seaice_alb = ', iflag_seaice_alb 2715 2720 WRITE(lunout,*) ' iflag_leads = ', iflag_leads 2721 WRITE(lunout,*) ' sic_hice_fixed = ', sic_hice_fixed 2716 2722 WRITE(lunout,*) ' sice_cond = ', sice_cond 2717 2723 WRITE(lunout,*) ' sisno_cond = ', sisno_cond … … 2730 2736 WRITE(lunout,*) ' fseaN = ', fseaN 2731 2737 WRITE(lunout,*) ' fseaS = ', fseaS 2732 !GG 2738 2733 2739 WRITE(lunout,*) ' n2o_cycle_cpl = ', n2o_cycle_cpl 2734 2740 WRITE(lunout,*) ' ndp_cycle_cpl = ', ndp_cycle_cpl -
LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90
r6053 r6070 4 4 MODULE ocean_forced_mod 5 5 ! 6 ! This module is used for both the sub-surfaces ocean and sea-ice for the case of a 6 ! This module is used for both the sub-surfaces ocean and sea-ice for the case of a 7 7 ! forced ocean, "ocean=force". 8 8 ! 9 IMPLICIT NONE9 IMPLICIT NONE 10 10 11 11 CONTAINS … … 13 13 !**************************************************************************************** 14 14 ! 15 SUBROUTINE ocean_forced_noice( &16 17 18 19 20 21 22 radsol, snow, agesno, &23 24 25 dthetadz300,pctsrf,Ampl &26 #ifdef ISO 27 ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &28 xtsnow,xtevap,h1 &29 #endif 30 15 SUBROUTINE ocean_forced_noice( & 16 itime, dtime, jour, knon, knindex, & 17 p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, & 18 temp_air, spechum, & 19 AcoefH, AcoefQ, BcoefH, BcoefQ, & 20 AcoefU, AcoefV, BcoefU, BcoefV, & 21 ps, u1, v1, gustiness, tsurf_in, & 22 radsol, snow, agesno, & 23 qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & 24 tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa, & 25 dthetadz300, pctsrf, Ampl & 26 #ifdef ISO 27 , xtprecip_rain, xtprecip_snow, xtspechum, Roce, rlat, & 28 xtsnow, xtevap, h1 & 29 #endif 30 ) 31 31 !$gpum horizontal knon klon 32 32 33 ! 34 ! This subroutine treats the "open ocean", all grid points that are not entierly covered 35 ! by ice. 36 ! The routine receives data from climatologie file limit.nc and does some calculations at the 37 ! surface. 38 ! 39 USE dimphy 40 USE calcul_fluxs_mod 41 USE limit_read_mod 42 USE mod_grid_phy_lmdz 43 USE indice_sol_mod 44 USE surface_data, ONLY : iflag_leads 45 USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o 46 use config_ocean_skin_m, only: activate_ocean_skin 47 USE calbeta_mod, ONLY : calbeta 48 49 #ifdef ISO 50 USE infotrac_phy, ONLY: ntiso,niso 51 USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall 33 !=========================================================================================== 34 ! This subroutine treats the "open ocean" tile in forced mode that is, when LMDZ is not 35 ! coupled to the slab or to NEMO. 36 ! The routine receives data from boundary-file limit.nc and does some calculations at the 37 ! surface. 38 !=========================================================================================== 39 40 USE dimphy 41 USE calcul_fluxs_mod 42 USE limit_read_mod 43 USE mod_grid_phy_lmdz 44 USE indice_sol_mod 45 USE surface_data, ONLY: iflag_leads 46 USE phys_output_var_mod, ONLY: sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o 47 use config_ocean_skin_m, only: activate_ocean_skin 48 USE calbeta_mod, ONLY: calbeta 49 50 #ifdef ISO 51 USE infotrac_phy, ONLY: ntiso, niso 52 USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall 52 53 #ifdef ISOVERIF 53 USE isotopes_mod, ONLY: iso_eau,ridicule54 !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix55 USE isotopes_verif_mod56 #endif 57 #endif 58 USE flux_arp_mod_h59 60 USE yomcst_mod_h54 USE isotopes_mod, ONLY: iso_eau, ridicule 55 !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix 56 USE isotopes_verif_mod 57 #endif 58 #endif 59 USE flux_arp_mod_h 60 USE clesphys_mod_h 61 USE yomcst_mod_h 61 62 62 63 ! Input arguments 63 64 !**************************************************************************************** 64 INTEGER, INTENT(IN) :: itime, jour, knon 65 INTEGER, DIMENSION(knon), INTENT(IN) :: knindex 66 REAL, INTENT(IN) :: dtime 67 REAL, DIMENSION(knon), INTENT(IN) :: p1lay 68 REAL, DIMENSION(knon), INTENT(IN) :: cdragh, cdragq, cdragm 69 REAL, DIMENSION(knon), INTENT(IN) :: precip_rain, precip_snow 70 REAL, DIMENSION(knon), INTENT(IN) :: temp_air, spechum 71 REAL, DIMENSION(knon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ 72 REAL, DIMENSION(knon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV 73 REAL, DIMENSION(knon), INTENT(IN) :: ps 74 REAL, DIMENSION(knon), INTENT(IN) :: u1, v1, gustiness 75 REAL, DIMENSION(knon), INTENT(IN) :: tsurf_in 76 real, intent(in):: rhoa(knon) ! (knon) density of moist air (kg / m3) 77 !GG 78 REAL, DIMENSION(knon), INTENT(IN) :: dthetadz300 79 REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf 80 ! 81 82 #ifdef ISO 83 REAL, DIMENSION(ntiso,knon), INTENT(IN) :: xtprecip_rain, xtprecip_snow 84 REAL, DIMENSION(ntiso,knon), INTENT(IN) :: xtspechum 85 REAL, DIMENSION(klon), INTENT(IN) :: rlat 65 INTEGER, INTENT(IN) :: itime, jour, knon 66 INTEGER, DIMENSION(knon), INTENT(IN) :: knindex 67 REAL, INTENT(IN) :: dtime 68 REAL, DIMENSION(knon), INTENT(IN) :: p1lay 69 REAL, DIMENSION(knon), INTENT(IN) :: cdragh, cdragq, cdragm 70 REAL, DIMENSION(knon), INTENT(IN) :: precip_rain, precip_snow 71 REAL, DIMENSION(knon), INTENT(IN) :: temp_air, spechum 72 REAL, DIMENSION(knon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ 73 REAL, DIMENSION(knon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV 74 REAL, DIMENSION(knon), INTENT(IN) :: ps 75 REAL, DIMENSION(knon), INTENT(IN) :: u1, v1, gustiness 76 REAL, DIMENSION(knon), INTENT(IN) :: tsurf_in 77 REAL, INTENT(IN) :: rhoa(knon) ! (knon) density of moist air (kg / m3) 78 REAL, DIMENSION(knon), INTENT(IN) :: dthetadz300 79 REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf 80 81 #ifdef ISO 82 REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtprecip_rain, xtprecip_snow 83 REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtspechum 84 REAL, DIMENSION(klon), INTENT(IN) :: rlat 86 85 #endif 87 86 88 87 ! In/Output arguments 89 88 !**************************************************************************************** 90 REAL, DIMENSION(knon), INTENT(INOUT) :: radsol91 REAL, DIMENSION(knon), INTENT(INOUT) :: snow92 REAL, DIMENSION(knon), INTENT(INOUT) :: agesno !? put to 0 in ocean93 #ifdef ISO 94 REAL, DIMENSION(niso,knon), INTENT(IN) :: xtsnow95 REAL, DIMENSION(niso,knon), INTENT(INOUT):: Roce89 REAL, DIMENSION(knon), INTENT(INOUT) :: radsol 90 REAL, DIMENSION(knon), INTENT(INOUT) :: snow 91 REAL, DIMENSION(knon), INTENT(INOUT) :: agesno !? put to 0 in ocean 92 #ifdef ISO 93 REAL, DIMENSION(niso, knon), INTENT(IN) :: xtsnow 94 REAL, DIMENSION(niso, knon), INTENT(INOUT):: Roce 96 95 #endif 97 96 98 97 ! Output arguments 99 98 !**************************************************************************************** 100 REAL, DIMENSION(knon), INTENT(OUT) :: qsurf 101 REAL, DIMENSION(knon), INTENT(OUT) :: evap, fluxsens, fluxlat 102 REAL, DIMENSION(knon), INTENT(OUT) :: flux_u1, flux_v1 103 REAL, DIMENSION(knon), INTENT(OUT) :: tsurf_new 104 REAL, DIMENSION(knon), INTENT(OUT) :: dflux_s, dflux_l 105 REAL, INTENT(out):: sens_prec_liq(:) ! (knon) 106 !GG 107 REAL, DIMENSION(knon), INTENT(OUT) :: Ampl 108 ! 109 110 #ifdef ISO 111 REAL, DIMENSION(ntiso,knon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux 112 REAL, DIMENSION(knon), INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation 99 REAL, DIMENSION(knon), INTENT(OUT) :: qsurf 100 REAL, DIMENSION(knon), INTENT(OUT) :: evap, fluxsens, fluxlat 101 REAL, DIMENSION(knon), INTENT(OUT) :: flux_u1, flux_v1 102 REAL, DIMENSION(knon), INTENT(OUT) :: tsurf_new 103 REAL, DIMENSION(knon), INTENT(OUT) :: dflux_s, dflux_l 104 REAL, INTENT(out):: sens_prec_liq(:) ! (knon) 105 REAL, DIMENSION(knon), INTENT(OUT) :: Ampl 106 107 #ifdef ISO 108 REAL, DIMENSION(ntiso, knon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux 109 REAL, DIMENSION(knon), INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation 113 110 #endif 114 111 115 112 ! Local variables 116 113 !**************************************************************************************** 117 INTEGER :: i, j 118 REAL, DIMENSION(knon) :: cal, beta, dif_grnd 119 REAL, DIMENSION(knon) :: alb_neig, tsurf_lim, zx_sl 120 REAL, DIMENSION(knon) :: u0, v0 121 REAL, DIMENSION(knon) :: u1_lay, v1_lay 122 LOGICAL :: check=.FALSE. 123 REAL, DIMENSION(knon) :: sens_prec_sol 124 REAL, DIMENSION(knon) :: lat_prec_liq, lat_prec_sol 125 ! GG 126 REAL, DIMENSION(knon) :: l_CBL, sicfra 127 ! 128 #ifdef ISO 129 REAL, PARAMETER :: t_coup = 273.15 130 #endif 131 114 INTEGER :: i, j 115 REAL, DIMENSION(knon) :: cal, beta, dif_grnd 116 REAL, DIMENSION(knon) :: alb_neig, tsurf_lim, zx_sl 117 REAL, DIMENSION(knon) :: u0, v0 118 REAL, DIMENSION(knon) :: u1_lay, v1_lay 119 LOGICAL :: check = .FALSE. 120 REAL, DIMENSION(knon) :: sens_prec_sol 121 REAL, DIMENSION(knon) :: lat_prec_liq, lat_prec_sol 122 REAL, DIMENSION(knon) :: l_CBL, sicfra 123 #ifdef ISO 124 REAL, PARAMETER :: t_coup = 273.15 125 #endif 132 126 133 127 !**************************************************************************************** 134 128 ! Start calculation 135 129 !**************************************************************************************** 136 IF (check) WRITE(*,*)' Entering ocean_forced_noice'130 IF (check) WRITE (*, *) ' Entering ocean_forced_noice' 137 131 138 132 #ifdef ISO 139 133 #ifdef ISOVERIF 140 DO i = 1, knon141 IF (iso_eau > 0) THEN142 CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), &143 & spechum(i),'ocean_forced_mod 111', &144 & errmax,errmaxrel)145 CALL iso_verif_egalite_choix(snow(i), &146 & xtsnow(iso_eau,i),'ocean_forced_mod 117', &147 & errmax,errmaxrel)148 ENDIF !IF (iso_eau > 0) THEN149 ENDDO !DO i=1,knon150 #endif 151 #endif 152 153 !**************************************************************************************** 154 ! 1) 134 DO i = 1, knon 135 IF (iso_eau > 0) THEN 136 CALL iso_verif_egalite_choix(xtspechum(iso_eau, i), & 137 & spechum(i), 'ocean_forced_mod 111', & 138 & errmax, errmaxrel) 139 CALL iso_verif_egalite_choix(snow(i), & 140 & xtsnow(iso_eau, i), 'ocean_forced_mod 117', & 141 & errmax, errmaxrel) 142 END IF !IF (iso_eau > 0) THEN 143 END DO !DO i=1,knon 144 #endif 145 #endif 146 147 !**************************************************************************************** 148 ! 1) 155 149 ! Read sea-surface temperature from file limit.nc 156 150 ! 157 151 !**************************************************************************************** 158 !--sb: 159 !!jyg if (knon.eq.1) then ! single-column model 160 if (klon_glo.eq.1) then ! single-column model 161 ! EV: now surface Tin flux_arp.h 162 !CALL read_tsurf1d(knon,tsurf_lim) ! new 163 DO i = 1, knon 164 tsurf_lim(i) = tg 165 ENDDO 166 167 else ! GCM 168 CALL limit_read_sst(knon,knindex,tsurf_lim & 169 #ifdef ISO 170 & ,Roce,rlat & 171 #endif 172 & ) 173 endif ! knon 174 !sb-- 152 IF (klon_glo .eq. 1) THEN ! single-column model 153 DO i = 1, knon 154 tsurf_lim(i) = tg 155 END DO 156 157 ELSE ! GCM 158 CALL limit_read_sst(knon, knindex, tsurf_lim & 159 #ifdef ISO 160 , Roce, rlat & 161 #endif 162 ) 163 END IF ! knon 175 164 176 165 !**************************************************************************************** … … 180 169 !**************************************************************************************** 181 170 ! Set some variables for calcul_fluxs 182 !cal = 0. 183 !beta = 1. 184 !dif_grnd = 0. 185 186 187 ! EV: use calbeta to calculate beta 188 ! Need to initialize qsurf for calbeta but it is not modified by this routine 189 qsurf(:)=0. 190 CALL calbeta(dtime, is_oce, knon, snow, qsurf, beta, cal, dif_grnd) 191 192 193 alb_neig(:) = 0. 194 agesno(:) = 0. 195 lat_prec_liq = 0.; lat_prec_sol = 0. 171 172 ! Use calbeta to calculate beta 173 ! Need to initialize qsurf for calbeta but it is not modified by this routine 174 qsurf(:) = 0. 175 CALL calbeta(dtime, is_oce, knon, snow, qsurf, beta, cal, dif_grnd) 176 177 alb_neig(:) = 0. 178 agesno(:) = 0. 179 lat_prec_liq = 0.; lat_prec_sol = 0. 196 180 197 181 ! Suppose zero surface speed 198 u0(:)=0.0 199 v0(:)=0.0 200 u1_lay(:) = u1(:) - u0(:) 201 v1_lay(:) = v1(:) - v0(:) 202 203 ! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf 204 IF (activate_ocean_skin == 2) THEN 205 CALL calcul_fluxs(knon, is_oce, dtime, & 206 tsurf_in, p1lay, cal, & 207 beta, cdragh, cdragq, ps, & 208 precip_rain, precip_snow, snow, qsurf, & 209 radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & 210 f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & 211 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & 212 sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) 182 u0(:) = 0.0 183 v0(:) = 0.0 184 u1_lay(:) = u1(:) - u0(:) 185 v1_lay(:) = v1(:) - v0(:) 186 187 ! Computation of tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf 188 ! EV: does that make sense to have activate_ocean_skin in ocean_forced? 189 ! in anycase; cal=dif=0 for ocean surfaces, so tsurf_new=tsurf_in. 190 ! the two options below should give the same result in principle. 191 192 IF (activate_ocean_skin == 2) THEN 193 CALL calcul_fluxs(knon, is_oce, dtime, & 194 tsurf_in, p1lay, cal, & 195 beta, cdragh, cdragq, ps, & 196 precip_rain, precip_snow, snow, qsurf, & 197 radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & 198 f_qsat_oce, AcoefH, AcoefQ, BcoefH, BcoefQ, & 199 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & 200 sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) 213 201 tsurf_new = tsurf_lim 214 ELSE 215 CALL calcul_fluxs(knon, is_oce, dtime, & 216 tsurf_lim, p1lay, cal, & 217 beta, cdragh, cdragq, ps, & 218 precip_rain, precip_snow, snow, qsurf, & 219 radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & 220 f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & 221 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & 222 sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) 223 ENDIF 224 225 do j = 1, knon 226 i = knindex(j) 227 sens_prec_liq_o(i,1) = sens_prec_liq(j) 228 sens_prec_sol_o(i,1) = sens_prec_sol(j) 229 lat_prec_liq_o(i,1) = lat_prec_liq(j) 230 lat_prec_sol_o(i,1) = lat_prec_sol(j) 231 enddo 232 233 !GG 234 235 if (iflag_leads == 1) then 202 ELSE 203 CALL calcul_fluxs(knon, is_oce, dtime, & 204 tsurf_lim, p1lay, cal, & 205 beta, cdragh, cdragq, ps, & 206 precip_rain, precip_snow, snow, qsurf, & 207 radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & 208 f_qsat_oce, AcoefH, AcoefQ, BcoefH, BcoefQ, & 209 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & 210 sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) 211 END IF 212 213 do j = 1, knon 214 i = knindex(j) 215 sens_prec_liq_o(i, 1) = sens_prec_liq(j) 216 sens_prec_sol_o(i, 1) = sens_prec_sol(j) 217 lat_prec_liq_o(i, 1) = lat_prec_liq(j) 218 lat_prec_sol_o(i, 1) = lat_prec_sol(j) 219 end do 220 221 IF (iflag_leads .GE. 1) THEN 236 222 !ym hyper faux, melange de champ compresse et non compresse, je rétabli... 237 223 !ym l_CBL = -52381.*dthetadz300 + 3008.1 … … 244 230 !ym fluxsens=Ampl*fluxsens 245 231 !ym dflux_s=Ampl*dflux_s 246 247 do j = 1, knon 248 i = knindex(j) 249 l_CBL(j) = -52381.*dthetadz300(j) + 3008.1 250 Ampl(j) = 6.012e-08*l_CBL(j)**2 - 4.036e-04*l_CBL(j) + 1.4979 251 IF (Ampl(j)>1.2) Ampl(j)=1.2 252 sicfra(j)=pctsrf(i,is_sic)/(1.-pctsrf(i,is_lic)-pctsrf(i,is_ter)) 253 IF (pctsrf(i,is_sic)+pctsrf(i,is_oce)<EPSFRA) sicfra(j)=0. 254 IF (sicfra(j)<0.7) Ampl(j)=1. 255 IF (sicfra(j)>0.7 .and. sicfra(j)<0.9) Ampl(j)=((sicfra(j)-0.7)/0.2)*Ampl(j)+((0.9-sicfra(j))/0.2) 256 fluxsens(j)=Ampl(j)*fluxsens(j) 257 dflux_s(j)=Ampl(j)*dflux_s(j) 258 enddo 259 endif 260 261 262 ! - Flux calculation at first modele level for U and V 263 CALL calcul_flux_wind(knon, dtime, & 264 u0, v0, u1, v1, gustiness, cdragm, & 265 AcoefU, AcoefV, BcoefU, BcoefV, & 266 p1lay, temp_air, & 267 flux_u1, flux_v1) 268 269 #ifdef ISO 270 CALL calcul_iso_surf_oce_vectall(knon, knon,t_coup, & 271 & ps,tsurf_new,spechum,u1_lay, v1_lay, xtspechum, & 272 & evap, Roce,xtevap,h1 & 232 233 ! Amplification of surface fluxes over leads regions by G. Gastineau 234 ! Ref: Tian Tian et al. 2025, The Cryosphere 235 ! Isau I. 2007, JGR 236 ! Davy and Gao 2019, zenodo 237 ! Only sensible heat flux is accounted for if iflag_leads=1, evaporative and latent 238 ! heat fluxes are also considered if >1. 239 240 DO j = 1, knon 241 i = knindex(j) 242 l_CBL(j) = -52381.*dthetadz300(j) + 3008.1 243 Ampl(j) = 6.012e-08*l_CBL(j)**2 - 4.036e-04*l_CBL(j) + 1.4979 244 IF (Ampl(j) > 1.2) Ampl(j) = 1.2 245 sicfra(j) = pctsrf(i, is_sic)/(1.-pctsrf(i, is_lic) - pctsrf(i, is_ter)) 246 IF (pctsrf(i, is_sic) + pctsrf(i, is_oce) < EPSFRA) sicfra(j) = 0. 247 IF (sicfra(j) < 0.7) Ampl(j) = 1. 248 IF (sicfra(j) > 0.7 .and. sicfra(j) < 0.9) Ampl(j) = ((sicfra(j) - 0.7)/0.2)*Ampl(j) + ((0.9 - sicfra(j))/0.2) 249 fluxsens(j) = Ampl(j)*fluxsens(j) 250 dflux_s(j) = Ampl(j)*dflux_s(j) 251 END DO 252 253 IF (iflag_leads .GE. 2) THEN 254 fluxlat(j) = Ampl(j)*fluxlat(j) 255 dflux_l(j) = Ampl(j)*dflux_l(j) 256 evap(j) = Ampl(j)*evap(j) 257 END IF 258 259 END IF 260 261 ! - Flux calculation at first model level for U and V 262 CALL calcul_flux_wind(knon, dtime, & 263 u0, v0, u1, v1, gustiness, cdragm, & 264 AcoefU, AcoefV, BcoefU, BcoefV, & 265 p1lay, temp_air, & 266 flux_u1, flux_v1) 267 268 #ifdef ISO 269 CALL calcul_iso_surf_oce_vectall(knon, knon, t_coup, & 270 ps, tsurf_new, spechum, u1_lay, v1_lay, xtspechum, & 271 evap, Roce, xtevap, h1 & 273 272 #ifdef ISOTRAC 274 & ,knindex &275 #endif 276 &)277 #endif 273 , knindex & 274 #endif 275 ) 276 #endif 278 277 279 278 #ifdef ISO 280 279 #ifdef ISOVERIF 281 ! write(*,*) 'ocean_forced_mod 176: sortie de ocean_forced_noice' 282 IF (iso_eau > 0) THEN 283 DO i = 1, knon 284 CALL iso_verif_egalite_choix(snow(i), & 285 & xtsnow(iso_eau,i),'ocean_forced_mod 180', & 286 & errmax,errmaxrel) 287 ENDDO ! DO j=1,knon 288 ENDIF !IF (iso_eau > 0) THEN 289 #endif 290 #endif 291 292 END SUBROUTINE ocean_forced_noice 280 IF (iso_eau > 0) THEN 281 DO i = 1, knon 282 CALL iso_verif_egalite_choix(snow(i), & 283 xtsnow(iso_eau, i), 'ocean_forced_mod 180', & 284 errmax, errmaxrel) 285 END DO ! DO j=1,knon 286 END IF !IF (iso_eau > 0) THEN 287 #endif 288 #endif 289 290 END SUBROUTINE ocean_forced_noice 293 291 ! 294 292 !*************************************************************************************** 295 293 ! 296 SUBROUTINE ocean_forced_ice( &297 298 tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air,spechum, &299 300 301 302 303 304 305 306 307 308 #ifdef ISO 309 ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &310 xtsnow, xtsol,xtevap,Rland_ice &311 #endif 312 294 SUBROUTINE ocean_forced_ice( & 295 itime, dtime, jour, knon, knindex, & 296 tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, & 297 AcoefH, AcoefQ, BcoefH, BcoefQ, & 298 AcoefU, AcoefV, BcoefU, BcoefV, & 299 ps, u1, v1, gustiness, pctsrf, & 300 radsol, snow, qsol, agesno, tsoil, & 301 qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & 302 icesub, icemelt, & 303 tsurf_new, dflux_s, dflux_l, rhoa, swnet, hice, tice, bilg_cumul, & 304 fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, & 305 dtice_melt, dtice_snow2sic & 306 #ifdef ISO 307 , xtprecip_rain, xtprecip_snow, xtspechum, Roce, & 308 xtsnow, xtsol, xtevap, Rland_ice & 309 #endif 310 ) 313 311 !$gpum horizontal knon klon 314 312 ! 315 ! This subroutine treats the ocean where there is ice.316 ! Th e routine reads data from climatologie file and does flux calculations at the317 ! surface.318 ! 319 USE dimphy 320 USE geometry_mod, ONLY: longitude,latitude 321 USE calcul_fluxs_mod322 !GG USE surface_data, ONLY : calice, calsno 323 USE surface_data, ONLY : calice, calsno, iflag_seaice, iflag_seaice_alb, &324 sice_cond, sisno_cond, sisno_den, sisno_min, sithick_min, sisno_wfact, &325 amax_s,amax_n, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, &326 si_pen_frac, si_pen_ext, fseaN, fseaS, iflag_leads327 328 USE geometry_mod, ONLY: longitude,latitude,latitude_deg 329 !GG 330 USE limit_read_mod331 USE simplehydrol_mod, ONLY: simplehydrol332 USE indice_sol_mod333 USE albsno_mod, ONLY: albsno334 USE soil_mod, ONLY: soil335 USE calbeta_mod, ONLY: calbeta336 USE phys_output_var_mod, ONLY: sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o337 #ifdef ISO 338 USE infotrac_phy, ONLY: niso, ntiso339 USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall313 !=========================================================================================== 314 ! This subroutine treats the sea ice tile in when LMDZ is not coupled to the slab or NEMO 315 ! The routine reads data from boundary-condition file limit.nc and does flux calculations 316 ! at the surface. 317 !=========================================================================================== 318 319 USE dimphy 320 USE geometry_mod, ONLY: longitude, latitude 321 USE calcul_fluxs_mod 322 USE surface_data, ONLY: calice, calsno, iflag_seaice, iflag_seaice_alb, & 323 sice_cond, sisno_cond, sisno_den, sisno_min, sithick_min, sisno_wfact, & 324 amax_s, amax_n, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, & 325 si_pen_frac, si_pen_ext, fseaN, fseaS, iflag_leads, sic_hice_fixed 326 327 USE geometry_mod, ONLY: longitude, latitude, latitude_deg 328 USE limit_read_mod 329 USE simplehydrol_mod, ONLY: simplehydrol 330 USE indice_sol_mod 331 USE albsno_mod, ONLY: albsno 332 USE soil_mod, ONLY: soil 333 USE calbeta_mod, ONLY: calbeta 334 USE phys_output_var_mod, ONLY: sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o 335 #ifdef ISO 336 USE infotrac_phy, ONLY: niso, ntiso 337 USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall 340 338 #ifdef ISOVERIF 341 USE isotopes_mod, ONLY: iso_eau,ridicule 342 !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix 343 USE isotopes_verif_mod 344 #endif 345 #endif 346 USE flux_arp_mod_h 347 USE clesphys_mod_h 348 USE yomcst_mod_h 349 USE dimsoil_mod_h, ONLY: nsoilmx 350 351 ! INCLUDE "indicesol.h" 352 339 USE isotopes_mod, ONLY: iso_eau, ridicule 340 !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix 341 USE isotopes_verif_mod 342 #endif 343 #endif 344 USE flux_arp_mod_h 345 USE clesphys_mod_h 346 USE yomcst_mod_h 347 USE dimsoil_mod_h, ONLY: nsoilmx 353 348 354 349 ! Input arguments 355 350 !**************************************************************************************** 356 INTEGER, INTENT(IN) :: itime, jour, knon 357 INTEGER, DIMENSION(knon), INTENT(IN) :: knindex 358 REAL, INTENT(IN) :: dtime 359 REAL, DIMENSION(knon), INTENT(IN) :: tsurf_in 360 REAL, DIMENSION(knon), INTENT(IN) :: p1lay 361 REAL, DIMENSION(knon), INTENT(IN) :: cdragh, cdragm 362 REAL, DIMENSION(knon), INTENT(IN) :: precip_rain, precip_snow 363 REAL, DIMENSION(knon), INTENT(IN) :: temp_air, spechum 364 REAL, DIMENSION(knon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ 365 REAL, DIMENSION(knon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV 366 REAL, DIMENSION(knon), INTENT(IN) :: ps 367 REAL, DIMENSION(knon), INTENT(IN) :: u1, v1, gustiness 368 real, intent(in):: rhoa(knon) ! (knon) density of moist air (kg / m3) 369 !GG 370 REAL, DIMENSION(knon), INTENT(IN) :: swnet 371 REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf 372 !GG 373 #ifdef ISO 374 REAL, DIMENSION(ntiso,knon), INTENT(IN) :: xtprecip_rain, xtprecip_snow 375 REAL, DIMENSION(ntiso,knon), INTENT(IN) :: xtspechum 376 REAL, DIMENSION(niso,knon), INTENT(IN) :: Roce 377 REAL, DIMENSION(niso,knon), INTENT(IN) :: Rland_ice 351 INTEGER, INTENT(IN) :: itime, jour, knon 352 INTEGER, DIMENSION(knon), INTENT(IN) :: knindex 353 REAL, INTENT(IN) :: dtime 354 REAL, DIMENSION(knon), INTENT(IN) :: tsurf_in 355 REAL, DIMENSION(knon), INTENT(IN) :: p1lay 356 REAL, DIMENSION(knon), INTENT(IN) :: cdragh, cdragm 357 REAL, DIMENSION(knon), INTENT(IN) :: precip_rain, precip_snow 358 REAL, DIMENSION(knon), INTENT(IN) :: temp_air, spechum 359 REAL, DIMENSION(knon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ 360 REAL, DIMENSION(knon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV 361 REAL, DIMENSION(knon), INTENT(IN) :: ps 362 REAL, DIMENSION(knon), INTENT(IN) :: u1, v1, gustiness 363 real, intent(in):: rhoa(knon) ! (knon) density of moist air (kg / m3) 364 REAL, DIMENSION(knon), INTENT(IN) :: swnet 365 REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf 366 #ifdef ISO 367 REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtprecip_rain, xtprecip_snow 368 REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtspechum 369 REAL, DIMENSION(niso, knon), INTENT(IN) :: Roce 370 REAL, DIMENSION(niso, knon), INTENT(IN) :: Rland_ice 378 371 #endif 379 372 380 373 ! In/Output arguments 381 374 !**************************************************************************************** 382 REAL, DIMENSION(knon), INTENT(INOUT) :: radsol 383 REAL, DIMENSION(knon), INTENT(INOUT) :: snow, qsol 384 REAL, DIMENSION(knon), INTENT(INOUT) :: agesno 385 REAL, DIMENSION(knon, nsoilmx), INTENT(INOUT) :: tsoil 386 !GG 387 REAL, DIMENSION(klon), INTENT(INOUT) :: hice !ym WARNING uncompressed 388 REAL, DIMENSION(klon), INTENT(INOUT) :: tice !ym WARNING uncompressed 389 REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul !ym WARNING uncompressed 390 REAL, DIMENSION(klon), INTENT(INOUT) :: fcds !ym WARNING uncompressed 391 REAL, DIMENSION(klon), INTENT(INOUT) :: fcdi !ym WARNING uncompressed 392 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_growth !ym WARNING uncompressed 393 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_melt !ym WARNING uncompressed 394 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt !ym WARNING uncompressed 395 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_snow2sic !ym WARNING uncompressed 396 REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_melt !ym WARNING uncompressed 397 REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_snow2sic !ym WARNING uncompressed 398 !GG 399 #ifdef ISO 400 REAL, DIMENSION(niso,knon), INTENT(INOUT) :: xtsnow 401 REAL, DIMENSION(niso,knon), INTENT(IN) :: xtsol 375 REAL, DIMENSION(knon), INTENT(INOUT) :: radsol 376 REAL, DIMENSION(knon), INTENT(INOUT) :: snow, qsol 377 REAL, DIMENSION(knon), INTENT(INOUT) :: agesno 378 REAL, DIMENSION(knon, nsoilmx), INTENT(INOUT) :: tsoil 379 REAL, DIMENSION(klon), INTENT(INOUT) :: hice !ym WARNING uncompressed 380 REAL, DIMENSION(klon), INTENT(INOUT) :: tice !ym WARNING uncompressed 381 REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul !ym WARNING uncompressed 382 REAL, DIMENSION(klon), INTENT(INOUT) :: fcds !ym WARNING uncompressed 383 REAL, DIMENSION(klon), INTENT(INOUT) :: fcdi !ym WARNING uncompressed 384 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_growth !ym WARNING uncompressed 385 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_melt !ym WARNING uncompressed 386 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt !ym WARNING uncompressed 387 REAL, DIMENSION(klon), INTENT(INOUT) :: dh_snow2sic !ym WARNING uncompressed 388 REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_melt !ym WARNING uncompressed 389 REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_snow2sic !ym WARNING uncompressed 390 391 #ifdef ISO 392 REAL, DIMENSION(niso, knon), INTENT(INOUT) :: xtsnow 393 REAL, DIMENSION(niso, knon), INTENT(IN) :: xtsol 402 394 #endif 403 395 404 396 ! Output arguments 405 397 !**************************************************************************************** 406 REAL, DIMENSION(knon), INTENT(OUT) :: qsurf407 REAL, DIMENSION(knon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval408 REAL, DIMENSION(knon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval409 REAL, DIMENSION(knon), INTENT(OUT) :: evap, fluxsens, fluxlat410 REAL, DIMENSION(knon), INTENT(OUT) :: flux_u1, flux_v1411 REAL, DIMENSION(knon), INTENT(OUT) :: tsurf_new412 REAL, DIMENSION(knon), INTENT(OUT) :: dflux_s, dflux_l413 REAL, DIMENSION(knon), INTENT(OUT) :: icesub, icemelt414 #ifdef ISO 415 REAL, DIMENSION(ntiso,knon), INTENT(OUT) :: xtevap416 #endif 398 REAL, DIMENSION(knon), INTENT(OUT) :: qsurf 399 REAL, DIMENSION(knon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval 400 REAL, DIMENSION(knon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval 401 REAL, DIMENSION(knon), INTENT(OUT) :: evap, fluxsens, fluxlat 402 REAL, DIMENSION(knon), INTENT(OUT) :: flux_u1, flux_v1 403 REAL, DIMENSION(knon), INTENT(OUT) :: tsurf_new 404 REAL, DIMENSION(knon), INTENT(OUT) :: dflux_s, dflux_l 405 REAL, DIMENSION(knon), INTENT(OUT) :: icesub, icemelt 406 #ifdef ISO 407 REAL, DIMENSION(ntiso, knon), INTENT(OUT) :: xtevap 408 #endif 417 409 418 410 ! Local variables 419 411 !**************************************************************************************** 420 LOGICAL,PARAMETER :: check=.FALSE. 421 INTEGER :: i, j 422 REAL :: zfra 423 REAL, PARAMETER :: t_grnd=271.35 424 REAL, DIMENSION(knon) :: cal, beta, dif_grnd, capsol 425 REAL, DIMENSION(knon) :: alb_neig, tsurf_tmp 426 REAL, DIMENSION(knon) :: soilcap, soilflux 427 REAL, DIMENSION(knon) :: u0, v0 428 REAL, DIMENSION(knon) :: u1_lay, v1_lay 429 REAL, DIMENSION(knon) :: sens_prec_liq, sens_prec_sol 430 REAL, DIMENSION(knon) :: lat_prec_liq, lat_prec_sol 431 REAL, DIMENSION(knon) :: yhice ! compressed version of hice 432 !GG 433 INTEGER :: ki 434 INTEGER :: cpl_pas 435 REAL, DIMENSION(knon) :: bilg 436 REAL, DIMENSION(knon) :: fsic 437 REAL, DIMENSION(knon) :: f_bot 438 REAL, PARAMETER :: latent_ice = 334.0e3 439 REAL, PARAMETER :: rau_ice = 917.0 440 REAL, PARAMETER :: kice=2.2 441 REAL :: f_cond, f_swpen, f_cond_s, f_cond_i 442 REAL :: ustar, uscap, ustau 443 ! for snow/ice albedo: 444 REAL :: alb_snow, alb_ice, alb_pond 445 REAL :: frac_snow, frac_ice, frac_pond 446 REAL :: z1_i, z2_i, z1_s, zlog ! height parameters 447 ! for ice melt / freeze 448 REAL :: e_melt, snow_evap, h_test 449 ! dhsic, dfsic change in ice mass, fraction. 450 REAL :: dhsic, dfsic, frac_mf 451 REAL :: fsea, amax 452 REAL :: hice_i, tice_i, fsic_new 453 ! snow and ice physical characteristics: 454 REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp 455 REAL, PARAMETER :: t_melt=273.15 ! melting ice temp 456 REAL :: sno_den!=sisno_den !mean snow density, kg/m3 457 REAL, PARAMETER :: ice_den=917. ! ice density 458 REAL, PARAMETER :: sea_den=1025. ! sea water density 459 REAL :: ice_cond!=sice_cond*ice_den !conductivity of ice 460 REAL :: sno_cond!=sisno_cond*sno_den ! conductivity of snow 461 REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice 462 REAL, PARAMETER :: sea_cap=3995. ! specific heat capacity, water 463 REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice 464 465 ! control of snow and ice cover & freeze / melt (heights converted to kg/m2) 466 REAL :: snow_min!=sisno_min*sno_den !critical snow height 5 cm 467 REAL :: snow_wfact!=sisno_wfact ! max fraction of falling snow blown into ocean 468 REAL, PARAMETER :: ice_frac_min=0.005 469 REAL :: h_ice_min!=sithick_min ! min ice thickness 470 ! below ice_thin, priority is melt lateral / grow height 471 ! ice_thin is also height of new ice 472 REAL, PARAMETER :: h_ice_max=7 ! max ice height 473 ! Ice thickness parameter for lateral growth 474 REAL, PARAMETER :: h_ice_thick=1.5 475 REAL, PARAMETER :: h_ice_thin=0.15 412 LOGICAL, PARAMETER :: check = .FALSE. 413 INTEGER :: i, j 414 REAL :: zfra 415 REAL, PARAMETER :: t_grnd = 271.35 416 REAL, DIMENSION(knon) :: cal, beta, dif_grnd, capsol 417 REAL, DIMENSION(knon) :: alb_neig, tsurf_tmp 418 REAL, DIMENSION(knon) :: soilcap, soilflux 419 REAL, DIMENSION(knon) :: u0, v0 420 REAL, DIMENSION(knon) :: u1_lay, v1_lay 421 REAL, DIMENSION(knon) :: sens_prec_liq, sens_prec_sol 422 REAL, DIMENSION(knon) :: lat_prec_liq, lat_prec_sol 423 REAL, DIMENSION(knon) :: yhice ! compressed version of hice 424 INTEGER :: ki 425 INTEGER :: cpl_pas 426 REAL, DIMENSION(knon) :: bilg 427 REAL, DIMENSION(knon) :: fsic 428 REAL, DIMENSION(knon) :: f_bot 429 REAL, PARAMETER :: latent_ice = 334.0e3 430 REAL, PARAMETER :: rau_ice = 917.0 431 REAL, PARAMETER :: kice = 2.2 432 REAL :: f_cond, f_swpen, f_cond_s, f_cond_i 433 REAL :: ustar, uscap, ustau 434 ! for snow/ice albedo: 435 REAL :: alb_snow, alb_ice, alb_pond 436 REAL :: frac_snow, frac_ice, frac_pond 437 REAL :: z1_i, z2_i, z1_s, zlog ! height parameters 438 ! for ice melt / freeze 439 REAL :: e_melt, snow_evap, h_test 440 ! dhsic, dfsic change in ice mass, fraction. 441 REAL :: dhsic, dfsic, frac_mf 442 REAL :: fsea, amax 443 REAL :: hice_i, tice_i, fsic_new 444 ! snow and ice physical characteristics: 445 REAL, PARAMETER :: t_freeze = 271.35 ! freezing sea water temp 446 REAL, PARAMETER :: t_melt = 273.15 ! melting ice temp 447 REAL :: sno_den!=sisno_den !mean snow density, kg/m3 448 REAL, PARAMETER :: ice_den = 917. ! ice density 449 REAL, PARAMETER :: sea_den = 1025. ! sea water density 450 REAL :: ice_cond!=sice_cond*ice_den !conductivity of ice 451 REAL :: sno_cond!=sisno_cond*sno_den ! conductivity of snow 452 REAL, PARAMETER :: ice_cap = 2067. ! specific heat capacity, snow and ice 453 REAL, PARAMETER :: sea_cap = 3995. ! specific heat capacity, water 454 REAL, PARAMETER :: ice_lat = 334000. ! freeze /melt latent heat snow and ice 455 456 ! control of snow and ice cover & freeze / melt (heights converted to kg/m2) 457 REAL :: snow_min!=sisno_min*sno_den !critical snow height 5 cm 458 REAL :: snow_wfact!=sisno_wfact ! max fraction of falling snow blown into ocean 459 REAL, PARAMETER :: ice_frac_min = 0.005 460 REAL :: h_ice_min!=sithick_min ! min ice thickness 461 ! below ice_thin, priority is melt lateral / grow height 462 ! ice_thin is also height of new ice 463 REAL, PARAMETER :: h_ice_max = 7 ! max ice height 464 ! Ice thickness parameter for lateral growth 465 REAL, PARAMETER :: h_ice_thick = 1.5 466 REAL, PARAMETER :: h_ice_thin = 0.15 476 467 477 468 ! albedo and radiation parameters 478 INTEGER :: iflag_sic_albedo469 INTEGER :: iflag_sic_albedo 479 470 ! albedo old or NEMO 480 REAL :: alb_sno_dry!=rn_alb_sdry !dry snow albedo 481 REAL :: alb_sno_wet!=rn_alb_smlt !wet snow albedo 482 REAL :: alb_ice_dry!=rn_alb_idry !dry thick ice 483 REAL :: alb_ice_wet!=rn_alb_imlt !melting thick ice 484 ! new (Toyoda 2020) albedo 485 ! Values for snow / ice, dry / melting, visible / near IR 486 REAL, PARAMETER :: alb_sdry_vis=0.98 487 REAL, PARAMETER :: alb_smlt_vis=0.88 488 REAL, PARAMETER :: alb_sdry_nir=0.7 489 REAL, PARAMETER :: alb_smlt_nir=0.55 490 REAL, PARAMETER :: alb_idry_vis=0.78 491 REAL, PARAMETER :: alb_imlt_vis=0.705 492 REAL, PARAMETER :: alb_idry_nir=0.36 493 REAL, PARAMETER :: alb_imlt_nir=0.285 494 REAL, PARAMETER :: h_ice_alb=0.5*ice_den ! height for full ice albedo 495 REAL, PARAMETER :: h_sno_alb=0.02*300 ! height for control of snow fraction 496 497 REAL :: pen_frac !=si_pen_frac !fraction of shortwave penetrating into the 498 ! ice (not snow). Should be visible only, not NIR 499 REAL :: pen_ext !=si_pen_ext !extinction length of penetrating shortwave (m-1) 500 REAl :: lon(knon), lat(knon) ! for indexation 501 502 ! HF from ocean below ice 503 ! REAL, PARAMETER :: fseaN=2.0 ! NH 504 ! REAL, PARAMETER :: fseaS=4.0 ! SH 505 !GG 506 507 #ifdef ISO 508 REAL, PARAMETER :: t_coup = 273.15 509 REAL, DIMENSION(knon) :: fq_fonte_diag 510 REAL, DIMENSION(knon) :: fqfonte_diag 511 REAL, DIMENSION(knon) :: snow_evap_diag 512 REAL, DIMENSION(knon) :: fqcalving_diag 513 REAL, DIMENSION(knon) :: run_off_lic_diag 514 REAL :: coeff_rel_diag 515 REAL :: max_eau_sol_diag 516 REAL, DIMENSION(knon) :: runoff_diag 517 INTEGER IXT 518 REAL, DIMENSION(niso,knon) :: xtsnow_prec, xtsol_prec 519 REAL, DIMENSION(knon) :: snow_prec, qsol_prec 520 #endif 521 522 !**************************************************************************************** 523 ! Start calculation 524 !**************************************************************************************** 525 IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon 526 527 !**************************************************************************************** 528 ! 1) 529 ! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, flux_u1, flux_v1 530 ! dflux_s, dflux_l and qsurf 531 !**************************************************************************************** 532 533 tsurf_tmp(:) = tsurf_in(:) 534 535 !GG 536 IF (iflag_seaice==0) THEN ! Old LMDZ sea ice surface 537 !GG 538 ! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal 539 CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd) 540 541 542 IF (soil_model) THEN 543 ! update tsoil and calculate soilcap and soilflux 544 lon(1:knon) = longitude(knindex(1:knon)) 545 lat(1:knon) = latitude(knindex(1:knon)) 546 CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, & 547 & lon, lat, tsoil,soilcap, soilflux) 548 cal(1:knon) = RCPD / soilcap(1:knon) 549 radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) 550 dif_grnd = 1.0 / tau_gl 551 ELSE 552 dif_grnd = 1.0 / tau_gl 553 cal = RCPD * calice 554 WHERE (snow > 0.0) cal = RCPD * calsno 555 ENDIF 556 557 !GG 558 ELSEIF (iflag_seaice==2) THEN 559 560 sno_den=sisno_den !mean snow density, kg/m3 561 ice_cond=sice_cond*ice_den !conductivity of ice 562 sno_cond=sisno_cond*sno_den ! conductivity of snow 563 snow_min=sisno_min*sno_den !critical snow height 5 cm 564 snow_wfact=sisno_wfact ! max fraction of falling snow blown into ocean 565 h_ice_min=sithick_min ! min ice thickness 566 alb_sno_dry=rn_alb_sdry !dry snow albedo 567 alb_sno_wet=rn_alb_smlt !wet snow albedo 568 alb_ice_dry=rn_alb_idry !dry thick ice 569 alb_ice_wet=rn_alb_imlt !melting thick ice 570 pen_frac=si_pen_frac !fraction of shortwave penetrating into the 571 pen_ext=si_pen_ext !extinction length of penetrating shortwave (m-1) 572 573 bilg(:)=0. 574 dif_grnd(:)=0. 575 beta(:) = 1. 576 cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour 577 578 ! Surface, snow-ice and ice-ocean fluxes. 579 ! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd) 580 581 !write(*,*) 'radsol 1',radsol(1:100) 582 DO i=1,knon 583 ki=knindex(i) 584 fsic(i) = pctsrf(ki,is_sic) 585 IF (snow(i).GT.snow_min) THEN 586 ! 1 / snow-layer heat capacity 587 cal(i)=2.*RCPD/(snow(i)*ice_cap) 588 ! adjustment time-scale of conductive flux 589 dif_grnd(i) = cal(i) * sno_cond / snow(i) / RCPD 590 ! for conductive flux 591 f_cond_s = sno_cond * (tice(ki)-t_freeze) / snow(i) 592 radsol(i) = radsol(i)+f_cond_s 593 ! all shortwave flux absorbed 594 f_swpen=0. 595 ELSE ! bare ice. 596 f_cond_s = 0. 597 ! 1 / ice-layer heat capacity 598 cal(i) = 2.*RCPD/(hice(ki)*ice_den*ice_cap) 599 ! adjustment time-scale of conductive flux 600 dif_grnd(i) = cal(i) * ice_cond / (hice(ki)*ice_den) / RCPD 601 ! penetrative shortwave flux... 602 f_swpen=swnet(i)*pen_frac*exp(-pen_ext*hice(ki)) 603 radsol(i) = radsol(i)-f_swpen 604 ! GG no conductive flux in this case? 605 END IF 606 bilg(i)=f_swpen-f_cond_s 607 END DO 608 609 endif 610 611 ! beta = 1.0 612 lat_prec_liq = 0.; lat_prec_sol = 0. 613 614 ! Suppose zero surface speed 615 u0(:)=0.0 616 v0(:)=0.0 617 u1_lay(:) = u1(:) - u0(:) 618 v1_lay(:) = v1(:) - v0(:) 619 CALL calcul_fluxs(knon, is_sic, dtime, & 620 tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, & 621 precip_rain, precip_snow, snow, qsurf, & 622 radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & 623 f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & 624 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & 625 sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) 626 do j = 1, knon 627 i = knindex(j) 628 sens_prec_liq_o(i,2) = sens_prec_liq(j) 629 sens_prec_sol_o(i,2) = sens_prec_sol(j) 630 lat_prec_liq_o(i,2) = lat_prec_liq(j) 631 lat_prec_sol_o(i,2) = lat_prec_sol(j) 632 enddo 471 REAL :: alb_sno_dry!=rn_alb_sdry !dry snow albedo 472 REAL :: alb_sno_wet!=rn_alb_smlt !wet snow albedo 473 REAL :: alb_ice_dry!=rn_alb_idry !dry thick ice 474 REAL :: alb_ice_wet!=rn_alb_imlt !melting thick ice 475 ! Toyoda 2020' albedo 476 ! Values for snow / ice, dry / melting, visible / near IR 477 REAL, PARAMETER :: alb_sdry_vis = 0.98 478 REAL, PARAMETER :: alb_smlt_vis = 0.88 479 REAL, PARAMETER :: alb_sdry_nir = 0.7 480 REAL, PARAMETER :: alb_smlt_nir = 0.55 481 REAL, PARAMETER :: alb_idry_vis = 0.78 482 REAL, PARAMETER :: alb_imlt_vis = 0.705 483 REAL, PARAMETER :: alb_idry_nir = 0.36 484 REAL, PARAMETER :: alb_imlt_nir = 0.285 485 REAL, PARAMETER :: h_ice_alb = 0.5*ice_den ! height for full ice albedo 486 REAL, PARAMETER :: h_sno_alb = 0.02*300 ! height for control of snow fraction 487 488 REAL :: pen_frac !=si_pen_frac !fraction of shortwave penetrating into the 489 ! ice (not snow). Should be visible only, not NIR 490 REAL :: pen_ext !=si_pen_ext !extinction length of penetrating shortwave (m-1) 491 REAl :: lon(knon), lat(knon) ! for indexation 492 493 #ifdef ISO 494 REAL, PARAMETER :: t_coup = 273.15 495 REAL, DIMENSION(knon) :: fq_fonte_diag 496 REAL, DIMENSION(knon) :: fqfonte_diag 497 REAL, DIMENSION(knon) :: snow_evap_diag 498 REAL, DIMENSION(knon) :: fqcalving_diag 499 REAL, DIMENSION(knon) :: run_off_lic_diag 500 REAL :: coeff_rel_diag 501 REAL :: max_eau_sol_diag 502 REAL, DIMENSION(knon) :: runoff_diag 503 INTEGER IXT 504 REAL, DIMENSION(niso, knon) :: xtsnow_prec, xtsol_prec 505 REAL, DIMENSION(knon) :: snow_prec, qsol_prec 506 #endif 507 508 !**************************************************************************************** 509 ! Initialisation 510 !**************************************************************************************** 511 IF (check) WRITE (*, *) 'Entering surface_seaice, knon=', knon 512 513 iflag_sic_albedo = iflag_seaice_alb 514 515 IF (iflag_seaice .GE. 1) THEN 516 sno_den = sisno_den !mean snow density, kg/m3 517 ice_cond = sice_cond*ice_den !conductivity of ice 518 sno_cond = sisno_cond*sno_den ! conductivity of snow 519 snow_min = sisno_min*sno_den !critical snow height 5 cm 520 snow_wfact = sisno_wfact ! max fraction of falling snow blown into ocean 521 h_ice_min = sithick_min ! min ice thickness 522 alb_sno_dry = rn_alb_sdry !dry snow albedo 523 alb_sno_wet = rn_alb_smlt !wet snow albedo 524 alb_ice_dry = rn_alb_idry !dry thick ice 525 alb_ice_wet = rn_alb_imlt !melting thick ice 526 pen_frac = si_pen_frac !fraction of shortwave penetrating into the 527 pen_ext = si_pen_ext !extinction length of penetrating shortwave (m-1) 528 529 cpl_pas = NINT(86400./dtime*1.0) ! once a day 530 531 END IF 532 533 !*************************************************************************************************** 534 ! 1) Heat diffusion in the snow/sea-ice, surface turbulent flux and surface temperature calculation 535 ! if iflag_seaice=0, we use the old-style method considering a constant inertia for 536 ! sea ice (or uperlying snow) and eventually using the soil routine to compute heat diffusion 537 ! in sea ice 538 ! if iflag_seaice>1, we compute the sea-ice thermal properties using the scheme of the slab 539 ! from Liu, Risi, Codron et al. 2021, nature comm. 540 ! if iflag_seaice=1, we read the sea-ice thickness from limit.nc 541 ! if iflag_seaice=2, sea-ice thickness is computed using a prognostic evolution model. 542 ! see also PhD thesis by N. Michalezyk (LOCEAN) 543 !*************************************************************************************************** 544 545 tsurf_tmp(:) = tsurf_in(:) 546 547 IF (iflag_seaice == 0) THEN ! Old LMDZ sea ice surface 548 ! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal 549 CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd) 550 551 IF (soil_model) THEN 552 ! update tsoil and calculate soilcap and soilflux 553 lon(1:knon) = longitude(knindex(1:knon)) 554 lat(1:knon) = latitude(knindex(1:knon)) 555 CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, & 556 & lon, lat, tsoil, soilcap, soilflux) 557 cal(1:knon) = RCPD/soilcap(1:knon) 558 radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) 559 dif_grnd = 1.0/tau_gl 560 ELSE 561 dif_grnd = 1.0/tau_gl 562 cal = RCPD*calice 563 WHERE (snow > 0.0) cal = RCPD*calsno 564 END IF 565 566 ELSEIF (iflag_seaice .GE. 1) THEN 567 ! if iflag_seaice>=1, advanced computation of snow and sea-ice thermal properties. 568 ! They become dependent on sea-ice thickness hice. 569 bilg(:) = 0. 570 dif_grnd(:) = 0. 571 beta(:) = 1. 572 cpl_pas = NINT(86400./dtime*1.0) ! once a day 573 574 DO i = 1, knon 575 ki = knindex(i) 576 fsic(i) = pctsrf(ki, is_sic) 577 IF (snow(i) .GT. snow_min) THEN 578 ! 1 / snow-layer heat capacity 579 cal(i) = 2.*RCPD/(snow(i)*ice_cap) 580 ! adjustment time-scale of conductive flux 581 dif_grnd(i) = cal(i)*sno_cond/snow(i)/RCPD 582 ! for conductive flux 583 f_cond_s = sno_cond*(tice(ki) - t_freeze)/snow(i) 584 radsol(i) = radsol(i) + f_cond_s 585 ! all shortwave flux absorbed 586 f_swpen = 0. 587 ELSE ! bare ice. 588 f_cond_s = 0. 589 ! 1 / ice-layer heat capacity 590 cal(i) = 2.*RCPD/(hice(ki)*ice_den*ice_cap) 591 ! adjustment time-scale of conductive flux 592 dif_grnd(i) = cal(i)*ice_cond/(hice(ki)*ice_den)/RCPD 593 ! penetrative shortwave flux... 594 f_swpen = swnet(i)*pen_frac*exp(-pen_ext*hice(ki)) 595 radsol(i) = radsol(i) - f_swpen 596 ! GG no conductive flux in this case? 597 END IF 598 bilg(i) = f_swpen - f_cond_s 599 END DO 600 601 END IF 602 603 lat_prec_liq(:) = 0.0 604 lat_prec_sol(:) = 0.0 605 ! Suppose zero surface speed 606 u0(:) = 0.0 607 v0(:) = 0.0 608 u1_lay(:) = u1(:) - u0(:) 609 v1_lay(:) = v1(:) - v0(:) 610 CALL calcul_fluxs(knon, is_sic, dtime, & 611 tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, & 612 precip_rain, precip_snow, snow, qsurf, & 613 radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & 614 f_qsat_oce, AcoefH, AcoefQ, BcoefH, BcoefQ, & 615 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & 616 sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) 617 618 DO j = 1, knon 619 i = knindex(j) 620 sens_prec_liq_o(i, 2) = sens_prec_liq(j) 621 sens_prec_sol_o(i, 2) = sens_prec_sol(j) 622 lat_prec_liq_o(i, 2) = lat_prec_liq(j) 623 lat_prec_sol_o(i, 2) = lat_prec_sol(j) 624 END DO 633 625 634 626 ! - Flux calculation at first modele level for U and V 635 CALL calcul_flux_wind(knon, dtime, & 636 u0, v0, u1, v1, gustiness, cdragm, & 637 AcoefU, AcoefV, BcoefU, BcoefV, & 638 p1lay, temp_air, & 639 flux_u1, flux_v1) 640 641 !**************************************************************************************** 642 ! 2) 643 ! Calculations due to snow and runoff 644 ! 645 !**************************************************************************************** 646 !GG 647 if (iflag_seaice==0) then 648 !GG 649 #ifdef ISO 650 ! verif 627 CALL calcul_flux_wind(knon, dtime, & 628 u0, v0, u1, v1, gustiness, cdragm, & 629 AcoefU, AcoefV, BcoefU, BcoefV, & 630 p1lay, temp_air, & 631 flux_u1, flux_v1) 632 633 ! - If iflag_seaice >=1, we also compute ice temperature and heat flux in the ice 634 635 IF (iflag_seaice .GE. 1) THEN 636 DO i = 1, knon 637 ki = knindex(i) 638 639 ! ocean-ice heat flux 640 fsea = fseaS 641 amax = amax_s 642 if (latitude(ki) > 0) THEN 643 fsea = fseaN 644 amax = amax_n 645 END IF 646 647 IF (snow(i) .GT. snow_min) THEN 648 ! snow conductive flux after calcul_fluxs (pos up) 649 f_cond_s = sno_cond*(tice(ki) - tsurf_new(i))/snow(i) 650 ! 1 / heat capacity and conductive timescale 651 uscap = 2./ice_cap/(snow(i) + hice(ki)*ice_den) 652 ustau = uscap*ice_cond/(hice(ki)*ice_den) 653 ! update ice temp 654 tice(ki) = (tice(ki) + dtime*(ustau*t_freeze - uscap*f_cond_s)) & 655 /(1 + dtime*ustau) 656 ELSE ! bare ice 657 tice(ki) = tsurf_new(i) 658 END IF 659 ! ice conductive flux (pos up) 660 f_cond_i = ice_cond*(t_freeze - tice(ki))/(hice(ki)*ice_den) 661 f_bot(i) = fsea - f_cond_i 662 fcdi(ki) = f_cond_i - fsea 663 fcds(i) = f_cond_s 664 !bilg(i) = bilg(i)+f_cond_i 665 END DO 666 667 END IF 668 669 !**************************************************************************************** 670 ! 2) Treatment of snow and ice hydrology and computation of sea ice thickness 671 ! two different treatments if iflag_seaice = 0 (old-style) or >0 (slab style) 672 !**************************************************************************************** 673 674 IF (iflag_seaice .EQ. 0) THEN 675 676 ! if iflag_seaice = 0, we fix the sea-ice thickness to a constant value 677 ! this only serves for albedo calculation if iflag_sic_albedo >=1 678 DO i = 1, knon 679 ki = knindex(i) 680 hice(ki) = sic_hice_fixed 681 END DO 682 683 #ifdef ISO 651 684 #ifdef ISOVERIF 652 DO i = 1, knon 653 IF (iso_eau > 0) THEN 654 IF (snow(i) > ridicule) THEN 655 CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), & 656 & 'interfsurf 964',errmax,errmaxrel) 657 ENDIF !IF ((snow(i) > ridicule)) THEN 658 ENDIF !IF (iso_eau > 0) THEN 659 ENDDO !DO i=1,knon 660 #endif 661 ! end verif 662 663 DO i = 1, knon 664 snow_prec(i) = snow(i) 665 DO ixt = 1, niso 666 xtsnow_prec(ixt,i) = xtsnow(ixt,i) 667 ENDDO !DO ixt=1,niso 668 ! initialisation: 669 fq_fonte_diag(i) = 0.0 670 fqfonte_diag(i) = 0.0 671 snow_evap_diag(i)= 0.0 672 ENDDO !DO i=1,knon 673 #endif 674 675 676 CALL simplehydrol( knon, is_sic, knindex, dtime, & 677 tsurf_tmp, precip_rain, precip_snow, & 678 snow, qsol, tsurf_new, evap, icesub, icemelt & 679 #ifdef ISO 680 & ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag & 681 & ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag & 682 #endif 683 & ) 684 685 686 #ifdef ISO 687 ! isotopes: tout est externalisé 685 DO i = 1, knon 686 IF (iso_eau > 0) THEN 687 IF (snow(i) > ridicule) THEN 688 CALL iso_verif_egalite_choix(xtsnow(iso_eau, i), snow(i), & 689 & 'interfsurf 964', errmax, errmaxrel) 690 END IF !IF ((snow(i) > ridicule)) THEN 691 END IF !IF (iso_eau > 0) THEN 692 END DO !DO i=1,knon 693 #endif 694 695 DO i = 1, knon 696 snow_prec(i) = snow(i) 697 DO ixt = 1, niso 698 xtsnow_prec(ixt, i) = xtsnow(ixt, i) 699 END DO !DO ixt=1,niso 700 ! initialisation: 701 fq_fonte_diag(i) = 0.0 702 fqfonte_diag(i) = 0.0 703 snow_evap_diag(i) = 0.0 704 END DO !DO i=1,knon 705 #endif 706 707 CALL simplehydrol(knon, is_sic, knindex, dtime, & 708 tsurf_tmp, precip_rain, precip_snow, & 709 snow, qsol, tsurf_new, evap, icesub, icemelt & 710 #ifdef ISO 711 , fq_fonte_diag, fqfonte_diag, snow_evap_diag, fqcalving_diag & 712 , max_eau_sol_diag, runoff_diag, run_off_lic_diag, coeff_rel_diag & 713 #endif 714 ) 715 716 #ifdef ISO 717 CALL calcul_iso_surf_sic_vectall(knon, knon, & 718 evap, snow_evap_diag, Tsurf_new, Roce, snow, & 719 fq_fonte_diag, fqfonte_diag, dtime, t_coup, & 720 precip_snow, xtprecip_snow, xtprecip_rain, snow_prec, xtsnow_prec, & 721 xtspechum, spechum, ps, & 722 xtevap, xtsnow, fqcalving_diag, & 723 knindex, is_sic, run_off_lic_diag, coeff_rel_diag, Rland_ice) 724 #ifdef ISOVERIF 725 726 IF (iso_eau > 0) THEN 727 DO i = 1, knon 728 CALL iso_verif_egalite_choix(snow(i), & 729 xtsnow(iso_eau, i), 'ocean_forced_mod 396', & 730 errmax, errmaxrel) 731 END DO ! DO j=1,knon 732 END IF !IF (iso_eau > 0) then 733 #endif 688 734 !#ifdef ISOVERIF 689 ! write(*,*) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall' 690 ! write(*,*) 'klon,knon=',klon,knon 691 !#endif 692 CALL calcul_iso_surf_sic_vectall(knon,knon, & 693 & evap,snow_evap_diag,Tsurf_new,Roce,snow, & 694 & fq_fonte_diag,fqfonte_diag,dtime,t_coup, & 695 & precip_snow,xtprecip_snow,xtprecip_rain, snow_prec,xtsnow_prec, & 696 & xtspechum,spechum,ps, & 697 & xtevap,xtsnow,fqcalving_diag, & 698 & knindex,is_sic,run_off_lic_diag,coeff_rel_diag,Rland_ice & 699 & ) 700 #ifdef ISOVERIF 701 !write(*,*) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall' 702 IF (iso_eau > 0) THEN 703 DO i = 1, knon 704 CALL iso_verif_egalite_choix(snow(i), & 705 & xtsnow(iso_eau,i),'ocean_forced_mod 396', & 706 & errmax,errmaxrel) 707 ENDDO ! DO j=1,knon 708 ENDIF !IF (iso_eau > 0) then 709 #endif 710 !#ifdef ISOVERIF 711 #endif 735 #endif 712 736 !#ifdef ISO 713 714 ! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno) 715 ! 716 CALL albsno(knon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) 717 718 WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0. 719 720 alb1_new(:) = 0.0 721 DO i=1, knon 722 zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0))) 723 alb1_new(i) = alb_neig(i) * zfra + 0.6 * (1.0-zfra) 724 ENDDO 725 726 alb2_new(:) = alb1_new(:) 727 728 !GG 729 else 730 731 DO i=1,knon 732 ki=knindex(i) 733 734 ! ocean-ice heat flux 735 fsea=fseaS 736 amax=amax_s 737 if (latitude(ki)>0) THEN 738 fsea=fseaN 739 amax=amax_n 740 ENDIF 741 742 IF (snow(i).GT.snow_min) THEN 743 ! snow conductive flux after calcul_fluxs (pos up) 744 f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i) 745 ! 1 / heat capacity and conductive timescale 746 uscap = 2. / ice_cap / (snow(i)+hice(ki)*ice_den) 747 ustau = uscap * ice_cond / (hice(ki)*ice_den) 748 ! update ice temp 749 tice(ki) = (tice(ki) + dtime*(ustau*t_freeze - uscap*f_cond_s)) & 750 / (1 + dtime*ustau) 751 ELSE ! bare ice 752 tice(ki)=tsurf_new(i) 753 ENDIF 754 ! ice conductive flux (pos up) 755 f_cond_i = ice_cond * (t_freeze-tice(ki)) / (hice(ki)*ice_den) 756 f_bot(i) = fsea - f_cond_i 757 fcdi(ki) = f_cond_i - fsea 758 fcds(i) = f_cond_s 759 !bilg(i) = bilg(i)+f_cond_i 760 END DO 761 762 !**************************************************************************************** 763 ! 2) Update snow and ice surface : thickness 764 !**************************************************************************************** 765 766 IF (iflag_seaice==1) THEN 767 ! Read from limit 768 !ym totally wrong since "hice" is an uncompressed field (klon) and limit_read_hice return a compressed field (knon) 769 !ym CALL limit_read_hice(knon,knindex,hice) 770 CALL limit_read_hice(knon, knindex, yhice) 771 DO i=1,knon 772 ki=knindex(i) 773 hice(ki) = yhice(i) 774 ENDDO 775 ENDIF 776 ! Formula Krinner et al. 1997 : h = (0.2 + 3.8(f_min**2))(1 + 2(f- f_min)) 777 778 779 780 DO i=1,knon 781 ki=knindex(i) 782 !ym WARNING : fsic uninitilazed, take initialisation similar to iflag_seaice==2 783 fsic(i) = pctsrf(ki,is_sic) 784 IF (precip_snow(i) > 0.) THEN 785 snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(i))) 786 END IF 787 ! snow and ice sublimation 788 IF (evap(i) > 0.) THEN 789 snow_evap = MIN (snow(i) / dtime, evap(i)) 790 snow(i) = snow(i) - snow_evap * dtime 791 snow(i) = MAX(0.0, snow(i)) 792 IF (iflag_seaice==2) THEN 793 hice(ki) = MAX(0.0,hice(ki)-(evap(i)-snow_evap)*dtime/ice_den) 794 ENDIF 795 ENDIF 796 ! Melt / Freeze snow from above if Tsurf>0 797 IF (tsurf_new(i).GT.t_melt) THEN 798 ! energy available for melting snow (in kg of melted snow /m2) 799 e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. & 800 /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i)) 801 ! remove snow 802 tice_i=tice(ki) 803 IF (snow(i).GT.e_melt) THEN 804 snow(i)=snow(i)-e_melt 805 tsurf_new(i)=t_melt 806 ELSE ! all snow is melted 807 ! add remaining heat flux to ice 808 e_melt=e_melt-snow(i) 809 tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*hice(ki)*ice_den) 810 tsurf_new(i)=tice(ki) 811 END IF 812 dtice_melt(ki)=(tice(ki)-tice_i)/dtime 813 END IF 814 ! Bottom melt / grow 815 ! bottom freeze if bottom flux (cond + oce-ice) <0 816 IF (iflag_seaice==2) THEN 817 IF (f_bot(i).LT.0) THEN 818 ! larger fraction of bottom growth 819 frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thick) & 820 / (h_ice_max-h_ice_thick))) 821 ! quantity of new ice (formed at mean ice temp) 822 e_melt= -f_bot(i) * dtime * fsic(i) & 823 / (ice_lat+ice_cap/2.*(t_freeze-tice(ki))) 824 ! first increase height to h_thick 825 dhsic=MAX(0.,MIN(h_ice_thick-hice(ki),e_melt/(fsic(i)*ice_den))) 826 hice_i=hice(ki) 827 hice(ki)=dhsic+hice(ki) 828 e_melt=e_melt-fsic(i)*dhsic 829 IF (e_melt.GT.0.) THEN 830 ! frac_mf fraction used for lateral increase 831 dfsic=MIN(amax-fsic(i),e_melt*frac_mf/ (hice(ki)*ice_den) ) 832 ! No lateral growth -> forced ocean 833 !fsic(ki)=fsic(ki)+dfsic 834 e_melt=e_melt-dfsic*hice(ki)*ice_den 835 ! rest used to increase height 836 hice(ki)=MIN(h_ice_max,hice(ki)+e_melt/( fsic(i) * ice_den ) ) 837 END IF 838 dh_basal_growth(ki)=(hice(ki)-hice_i)/dtime 839 840 ! melt from below if bottom flux >0 841 ELSE 842 ! larger fraction of lateral melt from warm ocean 843 frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thin) & 844 / (h_ice_thick-h_ice_thin))) 845 ! bring ice to freezing and melt from below 846 ! quantity of melted ice 847 e_melt= f_bot(i) * dtime * fsic(i) & 848 / (ice_lat+ice_cap/2.*(tice(ki)-t_freeze)) 849 ! first decrease height to h_thick 850 hice_i=hice(ki) 851 dhsic=MAX(0.,MIN(hice(ki)-h_ice_thick,e_melt/(fsic(i)*ice_den))) 852 hice(ki)=hice(ki)-dhsic 853 e_melt=e_melt-fsic(i)*dhsic*ice_den 854 855 IF (e_melt.GT.0) THEN 856 ! frac_mf fraction used for height decrease 857 dhsic=MAX(0.,MIN(hice(ki)-h_ice_min,e_melt/ice_den*frac_mf/fsic(i))) 858 hice(ki)=hice(ki)-dhsic 859 e_melt=e_melt-fsic(i)*dhsic*ice_den 860 ! rest used to decrease fraction (up to 0!) 861 dfsic=MIN(fsic(i),e_melt/(hice(ki)*ice_den)) 862 ! Remaining heat not used if everything melted 863 e_melt=e_melt-dfsic*hice(ki)*ice_den 864 ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime 865 END IF 866 dh_basal_melt(ki)=(hice(ki)-hice_i)/dtime 737 738 ELSE 739 740 IF (iflag_seaice .EQ. 1) THEN 741 ! if iflag_seaice=1 we read the thickness of the sea ice in the limit.nc file 742 !ym totally wrong since "hice" is an uncompressed field (klon) and limit_read_hice return a compressed field (knon) 743 !ym CALL limit_read_hice(knon,knindex,hice) 744 CALL limit_read_hice(knon, knindex, yhice) 745 DO i = 1, knon 746 ki = knindex(i) 747 hice(ki) = yhice(i) 748 END DO 867 749 END IF 868 END IF 869 870 ! melt ice from above if Tice>0 871 tice_i=tice(ki) 872 IF (tice(ki).GT.t_melt) THEN 873 IF (iflag_seaice==2) THEN 874 ! quantity of ice melted (kg/m2) 875 e_melt=MAX(hice(ki)*ice_den*(tice(ki)-t_melt)*ice_cap/2. & 876 /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0) 877 ! melt from above, height only 878 hice_i=hice(ki) 879 dhsic=MIN(hice(ki)-h_ice_min,e_melt/ice_den) 880 dh_top_melt(i)=dhsic 881 e_melt=e_melt-dhsic 882 IF (e_melt.GT.0) THEN 883 ! lateral melt if ice too thin 884 dfsic=MAX(fsic(i)-ice_frac_min,e_melt/(h_ice_min*ice_den)*fsic(i)) 885 ! if all melted do nothing with remaining heat 886 e_melt=MAX(0.,e_melt*fsic(i)-dfsic*h_ice_min*ice_den) 887 ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime 888 END IF 889 hice(ki)=hice(ki)-dhsic 890 dh_top_melt(ki)=(hice(ki)-hice_i)/dtime 891 ! surface temperature at melting point 892 END IF 893 tice(ki)=t_melt 894 tsurf_new(i)=t_melt 895 END IF 896 dtice_melt(ki)=dtice_melt(ki)+tice(ki)-tice_i 897 898 ! convert snow to ice if below floating line 899 h_test=(hice(ki)*ice_den+snow(i))-hice(ki)*sea_den 900 IF ((h_test.GT.0.).AND.(hice(ki).GT.h_ice_min)) THEN !snow under water 901 ! extra snow converted to ice (with added frozen sea water) 902 IF (iflag_seaice==2) THEN 903 dhsic=h_test/(sea_den-ice_den+sno_den) 904 hice(ki)=hice(ki)+dhsic 905 ENDIF 906 snow(i)=snow(i)-dhsic*sno_den 907 ! available energy (freeze sea water + bring to tice) 908 e_melt=dhsic*ice_den*(1.-sno_den/ice_den)*(ice_lat+ & 909 ice_cap/2.*(t_freeze-tice(ki))) 910 ! update ice temperature 911 tice_i=tice(ki) 912 tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+hice(ki)*ice_den) 913 IF (iflag_seaice==2) THEN 914 dh_snow2sic(ki)=dhsic/dtime 915 END IF 916 dtice_snow2sic(ki)=(tice(ki)-tice_i)/dtime 917 END IF 918 END DO 919 920 !write(*,*) 'hice 2',hice(1:100) 921 !write(*,*) 'tice 2',tice(1:100) 922 923 iflag_sic_albedo=iflag_seaice_alb 924 925 !******************************************************************************* 926 ! 3) cumulate ice-ocean fluxes, update tslab, lateral grow 927 !***********************************************o******************************* 928 !cumul fluxes. 929 DO j=i,knon 930 ki=knindex(i) 931 bilg_cumul(ki)=bilg_cumul(ki) + bilg(j)/float(cpl_pas) 932 ENDDO 750 !NB another idea could be testing the formula by Krinner et al. 1997 : h = (0.2 + 3.8(f_min**2))(1 + 2(f- f_min)) 751 752 DO i = 1, knon 753 754 ki = knindex(i) 755 !ym WARNING : fsic uninitilazed, take initialisation similar to iflag_seaice==2 756 fsic(i) = pctsrf(ki, is_sic) 757 IF (precip_snow(i) > 0.) THEN 758 snow(i) = snow(i) + precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(i))) 759 END IF 760 ! snow and ice sublimation 761 IF (evap(i) > 0.) THEN 762 snow_evap = MIN(snow(i)/dtime, evap(i)) 763 snow(i) = snow(i) - snow_evap*dtime 764 snow(i) = MAX(0.0, snow(i)) 765 IF (iflag_seaice == 2) THEN 766 hice(ki) = MAX(0.0, hice(ki) - (evap(i) - snow_evap)*dtime/ice_den) 767 END IF 768 END IF 769 ! Melt / Freeze snow from above if Tsurf>0 770 IF (tsurf_new(i) .GT. t_melt) THEN 771 ! energy available for melting snow (in kg of melted snow /m2) 772 e_melt = MIN(MAX(snow(i)*(tsurf_new(i) - t_melt)*ice_cap/2. & 773 /(ice_lat + ice_cap/2.*(t_melt - tice(ki))), 0.0), snow(i)) 774 ! remove snow 775 tice_i = tice(ki) 776 IF (snow(i) .GT. e_melt) THEN 777 snow(i) = snow(i) - e_melt 778 tsurf_new(i) = t_melt 779 ELSE ! all snow is melted 780 ! add remaining heat flux to ice 781 e_melt = e_melt - snow(i) 782 tice(ki) = tice(ki) + e_melt*ice_lat*2./(ice_cap*hice(ki)*ice_den) 783 tsurf_new(i) = tice(ki) 784 END IF 785 dtice_melt(ki) = (tice(ki) - tice_i)/dtime 786 END IF 787 788 IF (iflag_seaice .EQ. 2) THEN 789 ! Interactive calculation of hice if iflag_seaice ==2 790 ! Bottom melt / grow 791 ! bottom freeze if bottom flux (cond + oce-ice) <0 792 793 IF (f_bot(i) .LT. 0) THEN 794 ! larger fraction of bottom growth 795 frac_mf = MIN(1., MAX(0., (hice(ki) - h_ice_thick) & 796 /(h_ice_max - h_ice_thick))) 797 ! quantity of new ice (formed at mean ice temp) 798 e_melt = -f_bot(i)*dtime*fsic(i) & 799 /(ice_lat + ice_cap/2.*(t_freeze - tice(ki))) 800 ! first increase height to h_thick 801 dhsic = MAX(0., MIN(h_ice_thick - hice(ki), e_melt/(fsic(i)*ice_den))) 802 hice_i = hice(ki) 803 hice(ki) = dhsic + hice(ki) 804 e_melt = e_melt - fsic(i)*dhsic 805 IF (e_melt .GT. 0.) THEN 806 ! frac_mf fraction used for lateral increase 807 dfsic = MIN(amax - fsic(i), e_melt*frac_mf/(hice(ki)*ice_den)) 808 ! No lateral growth -> forced ocean 809 !fsic(ki)=fsic(ki)+dfsic 810 e_melt = e_melt - dfsic*hice(ki)*ice_den 811 ! rest used to increase height 812 hice(ki) = MIN(h_ice_max, hice(ki) + e_melt/(fsic(i)*ice_den)) 813 END IF 814 dh_basal_growth(ki) = (hice(ki) - hice_i)/dtime 815 816 ! melt from below if bottom flux >0 817 ELSE 818 ! larger fraction of lateral melt from warm ocean 819 frac_mf = MIN(1., MAX(0., (hice(ki) - h_ice_thin) & 820 /(h_ice_thick - h_ice_thin))) 821 ! bring ice to freezing and melt from below 822 ! quantity of melted ice 823 e_melt = f_bot(i)*dtime*fsic(i) & 824 /(ice_lat + ice_cap/2.*(tice(ki) - t_freeze)) 825 ! first decrease height to h_thick 826 hice_i = hice(ki) 827 dhsic = MAX(0., MIN(hice(ki) - h_ice_thick, e_melt/(fsic(i)*ice_den))) 828 hice(ki) = hice(ki) - dhsic 829 e_melt = e_melt - fsic(i)*dhsic*ice_den 830 831 IF (e_melt .GT. 0) THEN 832 ! frac_mf fraction used for height decrease 833 dhsic = MAX(0., MIN(hice(ki) - h_ice_min, e_melt/ice_den*frac_mf/fsic(i))) 834 hice(ki) = hice(ki) - dhsic 835 e_melt = e_melt - fsic(i)*dhsic*ice_den 836 ! rest used to decrease fraction (up to 0!) 837 dfsic = MIN(fsic(i), e_melt/(hice(ki)*ice_den)) 838 ! Remaining heat not used if everything melted 839 e_melt = e_melt - dfsic*hice(ki)*ice_den 840 ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime 841 END IF 842 dh_basal_melt(ki) = (hice(ki) - hice_i)/dtime 843 END IF 844 END IF 845 846 ! melt ice from above if Tice>0 847 tice_i = tice(ki) 848 IF (tice(ki) .GT. t_melt) THEN 849 850 IF (iflag_seaice .EQ. 2) THEN 851 ! quantity of ice melted (kg/m2) 852 e_melt = MAX(hice(ki)*ice_den*(tice(ki) - t_melt)*ice_cap/2. & 853 /(ice_lat + ice_cap/2.*(t_melt - t_freeze)), 0.0) 854 ! melt from above, height only 855 hice_i = hice(ki) 856 dhsic = MIN(hice(ki) - h_ice_min, e_melt/ice_den) 857 dh_top_melt(i) = dhsic 858 e_melt = e_melt - dhsic 859 IF (e_melt .GT. 0) THEN 860 ! lateral melt if ice too thin 861 dfsic = MAX(fsic(i) - ice_frac_min, e_melt/(h_ice_min*ice_den)*fsic(i)) 862 ! if all melted do nothing with remaining heat 863 e_melt = MAX(0., e_melt*fsic(i) - dfsic*h_ice_min*ice_den) 864 ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime 865 END IF 866 hice(ki) = hice(ki) - dhsic 867 dh_top_melt(ki) = (hice(ki) - hice_i)/dtime 868 ! surface temperature at melting point 869 END IF 870 871 tice(ki) = t_melt 872 tsurf_new(i) = t_melt 873 END IF 874 dtice_melt(ki) = dtice_melt(ki) + tice(ki) - tice_i 875 876 ! convert snow to ice if below floating line 877 h_test = (hice(ki)*ice_den + snow(i)) - hice(ki)*sea_den 878 IF ((h_test .GT. 0.) .AND. (hice(ki) .GT. h_ice_min)) THEN !snow under water 879 ! extra snow converted to ice (with added frozen sea water) 880 IF (iflag_seaice .EQ. 2) THEN 881 dhsic = h_test/(sea_den - ice_den + sno_den) 882 hice(ki) = hice(ki) + dhsic 883 END IF 884 snow(i) = snow(i) - dhsic*sno_den 885 ! available energy (freeze sea water + bring to tice) 886 e_melt = dhsic*ice_den*(1.-sno_den/ice_den)*(ice_lat + & 887 ice_cap/2.*(t_freeze - tice(ki))) 888 ! update ice temperature 889 tice_i = tice(ki) 890 tice(ki) = tice(ki) + 2.*e_melt/ice_cap/(snow(i) + hice(ki)*ice_den) 891 IF (iflag_seaice .EQ. 2) THEN 892 dh_snow2sic(ki) = dhsic/dtime 893 END IF 894 dtice_snow2sic(ki) = (tice(ki) - tice_i)/dtime 895 END IF 896 897 END DO 898 899 ! cumulated ice-ocean fluxes, update tslab, lateral grow 900 DO j = i, knon 901 ki = knindex(i) 902 bilg_cumul(ki) = bilg_cumul(ki) + bilg(j)/float(cpl_pas) 903 END DO 933 904 934 905 !YM Pb for GPU port, done on an compressed field inside a compressed Kernel … … 938 909 !YM bilg_cumul(1:klon)=0. 939 910 !YM END IF ! coupling time 940 941 ! write(*,*) 'hice 3',hice(1:100) 942 ! write(*,*) 'tice 3',tice(1:100) 943 !tests ice fraction 911 !tests ice fraction 944 912 !YM Pb for GPU port : done on uncompressed field inside a compressed kernel 945 913 !YM update fsic only on compressed index … … 949 917 !YM END WHERE 950 918 951 DO j=i,knon 952 ki=knindex(i) 953 IF (fsic(i).LT.ice_frac_min ) THEN 954 tice(ki)=t_melt 955 hice(ki)=h_ice_min 956 ENDIF 957 ENDDO 958 959 !write(*,*) 'hice 4',hice(1:100) 960 !write(*,*) 'tice 4',tice(1:100) 961 962 endif 963 964 !**************************************************************************************** 965 ! 4) Compute sea-ice and snow albedo 966 !**************************************************************************************** 967 IF (iflag_seaice > 0) THEN 968 SELECT CASE (iflag_sic_albedo) 969 CASE(0) 970 ! old slab albedo : single value. age of snow, melt ponds. 971 DO i=1,knon 972 ki=knindex(i) 973 ! snow albedo: update snow age 974 IF (snow(i).GT.0.0001) THEN 975 agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)& 976 * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.) 977 ELSE 978 agesno(i)=0.0 979 END IF 980 ! snow albedo 981 alb_snow=alb_sno_wet+(alb_sno_dry-alb_sno_wet)*EXP(-agesno(i)/50.) 982 ! ice albedo (varies with ice tkickness and temp) 983 alb_ice=MAX(0.0,0.13*LOG(100.*hice(ki))+0.1) 984 !alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1) 985 IF (tice(ki).GT.t_freeze-0.01) THEN 986 alb_ice=MIN(alb_ice,alb_ice_wet) 987 ELSE 988 alb_ice=MIN(alb_ice,alb_ice_dry) 989 END IF 990 ! pond albedo 991 alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) 992 ! pond fraction 993 frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) 994 ! snow fraction 995 frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min)) 996 ! ice fraction 997 frac_ice=MAX(0.0,1.-frac_pond-frac_snow) 998 ! total albedo 999 alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond 1000 END DO 1001 alb2_new(:) = alb1_new(:) 1002 1003 CASE(1) 1004 ! New visible and IR albedos, dry / melting snow 1005 ! based on Toyoda et al, 2020 1006 DO i=1,knon 1007 ki=knindex(i) 1008 ! snow fraction 1009 frac_snow = snow(i) / (snow(i) + h_sno_alb) 1010 ! dependence of ice albedo with ice thickness 1011 frac_ice = MIN(1.,ATAN(4.*hice(ki)*ice_den) / ATAN(4.*h_ice_alb)) 1012 ! Total (for ice, min = 0.066 = alb_ocean) 1013 IF (tice(ki).GT.t_melt) THEN 1014 alb_ice = 0.066 + (alb_imlt_vis - 0.066)*frac_ice 1015 alb1_new(i)=alb_smlt_vis*frac_snow + alb_ice*(1.-frac_snow) 1016 alb_ice = 0.066 + (alb_imlt_nir - 0.066)*frac_ice 1017 alb2_new(i)=alb_smlt_nir*frac_snow + alb_ice*(1.-frac_snow) 1018 ELSEIF (tice(ki).GT.t_melt - 1.) THEN 1019 frac_pond = tice(ki) - t_freeze 1020 alb_snow = alb_smlt_vis*frac_pond + alb_sdry_vis*(1.-frac_pond) 1021 alb_ice = alb_imlt_vis*frac_pond + alb_idry_vis*(1.-frac_pond) 1022 alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice 1023 alb1_new(i)= alb_snow*frac_snow + alb_ice*(1.-frac_snow) 1024 alb_snow = alb_smlt_nir*frac_pond + alb_sdry_nir*(1.-frac_pond) 1025 alb_ice = alb_imlt_nir*frac_pond + alb_idry_nir*(1.-frac_pond) 1026 alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice 1027 alb2_new(i)= alb_snow*frac_snow + alb_ice*(1.-frac_snow) 1028 ELSE 1029 alb_ice = 0.066 + (alb_idry_vis - 0.066)*frac_ice 1030 alb1_new(i)=alb_sdry_vis*frac_snow + alb_ice*(1.-frac_snow) 1031 alb_ice = 0.066 + (alb_idry_nir - 0.066)*frac_ice 1032 alb2_new(i)=alb_sdry_nir*frac_snow + alb_ice*(1.-frac_snow) 1033 ENDIF 1034 END DO 1035 1036 CASE(2) 1037 ! LIM3 scheme. Uses clear sky / overcast value, with 50% clear sky 1038 z1_i = 1.5 * ice_den 1039 z2_i = 0.05 * ice_den 1040 zlog = 1. / (LOG(1.5) - LOG(0.05)) 1041 z1_s = 1. / (0.025 * sno_den) 1042 DO i=1,knon 1043 ki=knindex(i) 919 DO j = i, knon 920 ki = knindex(i) 921 IF (fsic(i) .LT. ice_frac_min) THEN 922 tice(ki) = t_melt 923 hice(ki) = h_ice_min 924 END IF 925 END DO 926 927 END IF ! if iflag_seaice .eq. 0 928 929 !**************************************************************************************** 930 ! 3) Compute sea-ice / snow albedo 931 ! there are currently 4 options 932 !**************************************************************************************** 933 934 SELECT CASE (iflag_sic_albedo) 935 CASE (0) 936 ! default option. Calculation of albedo at snow (alb_neig) and update the age of snow (agesno) 937 ! Note that where the is no snow, the albedo of bare sea ice is taken as the 938 ! value of aged snow (computation with agesno=0). 939 ! when looking at albsno, it is in fact a typical value of a desert (0.55)! 940 CALL albsno(knon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) 941 WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0. 942 alb1_new(:) = 0.0 943 DO i = 1, knon 944 zfra = MAX(0.0, MIN(1.0, snow(i)/(snow(i) + 10.0))) 945 alb1_new(i) = alb_neig(i)*zfra + 0.6*(1.0 - zfra) 946 END DO 947 alb2_new(:) = alb1_new(:) 948 949 CASE (1) 950 ! New visible and IR albedos, dry / melting snow 951 ! based on Toyoda et al, 2020, from the slab model 952 DO i = 1, knon 953 ki = knindex(i) 954 ! snow fraction 955 frac_snow = snow(i)/(snow(i) + h_sno_alb) 956 ! dependence of ice albedo with ice thickness 957 frac_ice = MIN(1., ATAN(4.*hice(ki)*ice_den)/ATAN(4.*h_ice_alb)) 958 ! Total (for ice, min = 0.066 = alb_ocean) 959 IF (tice(ki) .GT. t_melt) THEN 960 alb_ice = 0.066 + (alb_imlt_vis - 0.066)*frac_ice 961 alb1_new(i) = alb_smlt_vis*frac_snow + alb_ice*(1.-frac_snow) 962 alb_ice = 0.066 + (alb_imlt_nir - 0.066)*frac_ice 963 alb2_new(i) = alb_smlt_nir*frac_snow + alb_ice*(1.-frac_snow) 964 ELSEIF (tice(ki) .GT. t_melt - 1.) THEN 965 frac_pond = tice(ki) - t_freeze 966 alb_snow = alb_smlt_vis*frac_pond + alb_sdry_vis*(1.-frac_pond) 967 alb_ice = alb_imlt_vis*frac_pond + alb_idry_vis*(1.-frac_pond) 968 alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice 969 alb1_new(i) = alb_snow*frac_snow + alb_ice*(1.-frac_snow) 970 alb_snow = alb_smlt_nir*frac_pond + alb_sdry_nir*(1.-frac_pond) 971 alb_ice = alb_imlt_nir*frac_pond + alb_idry_nir*(1.-frac_pond) 972 alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice 973 alb2_new(i) = alb_snow*frac_snow + alb_ice*(1.-frac_snow) 974 ELSE 975 alb_ice = 0.066 + (alb_idry_vis - 0.066)*frac_ice 976 alb1_new(i) = alb_sdry_vis*frac_snow + alb_ice*(1.-frac_snow) 977 alb_ice = 0.066 + (alb_idry_nir - 0.066)*frac_ice 978 alb2_new(i) = alb_sdry_nir*frac_snow + alb_ice*(1.-frac_snow) 979 END IF 980 END DO 981 982 CASE (2) 983 ! LIM3 scheme. Uses clear sky / overcast value, assuming 50% clear sky 984 ! note that this is what is implemented in the coupled model (CMIP5, 6, 7) 985 z1_i = 1.5*ice_den 986 z2_i = 0.05*ice_den 987 zlog = 1./(LOG(1.5) - LOG(0.05)) 988 z1_s = 1./(0.025*sno_den) 989 DO i = 1, knon 990 ki = knindex(i) 1044 991 ! temperature above / below 0 1045 IF (tice(ki) .GE.t_melt) THEN992 IF (tice(ki) .GE. t_melt) THEN 1046 993 alb_ice = alb_ice_wet 1047 994 alb_snow = alb_sno_wet … … 1049 996 alb_ice = alb_ice_dry 1050 997 alb_snow = alb_sno_dry 1051 END IF998 END IF 1052 999 ! ice thickness 1053 IF (hice(ki)*ice_den .LT.z2_i) THEN1054 alb_ice = 0.066 + 0.114 * hice(ki)*ice_den /z2_i1055 ELSEIF (hice(ki)*ice_den .LT.z1_i) THEN1056 1057 * (LOG(z1_i) - LOG(hice(ki)*ice_den)) *zlog1058 END IF1000 IF (hice(ki)*ice_den .LT. z2_i) THEN 1001 alb_ice = 0.066 + 0.114*hice(ki)*ice_den/z2_i 1002 ELSEIF (hice(ki)*ice_den .LT. z1_i) THEN 1003 alb_ice = alb_ice + (0.18 - alb_ice) & 1004 *(LOG(z1_i) - LOG((hice(ki)+0.001)*ice_den))*zlog 1005 END IF 1059 1006 ! ice or snow depending on snow thickness 1060 alb_snow = alb_snow - (alb_snow - alb_ice) * EXP(- snow(i) *z1_s)1007 alb_snow = alb_snow - (alb_snow - alb_ice)*EXP(-snow(i)*z1_s) 1061 1008 ! Effect of clouds (polynomial fit with 50% clouds) 1062 alb1_new(i) = alb_snow - 0.5 * (-0.1010 *alb_snow*alb_snow &1063 + 0.1933*alb_snow - 0.0148)1009 alb1_new(i) = alb_snow - 0.5*(-0.1010*alb_snow*alb_snow & 1010 + 0.1933*alb_snow - 0.0148) 1064 1011 alb2_new(i) = alb1_new(i) 1065 END DO 1066 1067 CASE(3) 1068 CALL albsno(knon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) 1069 WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0. 1070 alb1_new(:) = 0.0 1071 DO i=1, knon 1072 zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0))) 1073 alb1_new(i) = alb_neig(i) * zfra + 0.6 * (1.0-zfra) 1074 ENDDO 1075 alb2_new(:) = alb1_new(:) 1076 1077 print*,'alb_neig=',alb_neig 1078 print*,'zfra=',zfra 1079 print*,'snow=',snow 1080 print*,'alb1_new=',alb1_new 1081 print*,'alb2_new=',alb2_new 1082 END SELECT 1083 END IF 1084 ! ------ End Albedo ---------- 1085 1086 !GG 1087 END SUBROUTINE ocean_forced_ice 1088 1089 SUBROUTINE ocean_forced_ice_reset_bilg_cumul(itime, dtime, bilg_cumul) 1090 USE dimphy, ONLY : klon 1091 USE surface_data, ONLY : iflag_seaice, type_ocean, version_ocean 1092 IMPLICIT NONE 1093 INTEGER, INTENT(IN) :: itime 1094 REAL, INTENT(IN) :: dtime 1095 REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul 1096 INTEGER :: cpl_pas 1097 1098 IF (iflag_seaice/=0) THEN 1099 cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour 1100 IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab 1101 IF ( .NOT. (type_ocean == 'couple' .OR. (type_ocean == 'slab'.AND.version_ocean=='sicINT'))) THEN 1102 bilg_cumul(1:klon)=0. 1103 ENDIF 1104 END IF ! coupling time 1105 ENDIF 1106 1107 END SUBROUTINE ocean_forced_ice_reset_bilg_cumul 1108 !************************************************************************ 1109 ! 1D case 1110 !************************************************************************ 1111 ! SUBROUTINE read_tsurf1d(knon,sst_out) 1112 ! 1113 ! This subroutine specifies the surface temperature to be used in 1D simulations 1114 ! 1115 ! USE dimphy, ONLY : klon 1116 ! 1117 ! INTEGER, INTENT(IN) :: knon ! nomber of points on compressed grid 1118 ! REAL, DIMENSION(klon), INTENT(OUT) :: sst_out ! tsurf used to force the single-column model 1119 ! 1120 ! INTEGER :: i 1121 ! COMMON defined in lmdz1d.F: 1122 ! real ts_cur 1123 ! common /sst_forcing/ts_cur 1124 ! 1125 ! DO i = 1, knon 1126 ! sst_out(i) = ts_cur 1127 ! ENDDO 1128 ! 1129 ! END SUBROUTINE read_tsurf1d 1130 ! 1131 ! 1132 !************************************************************************ 1012 END DO 1013 1014 CASE (3) 1015 ! old slab albedo : single value. age of snow, melt ponds. 1016 DO i = 1, knon 1017 ki = knindex(i) 1018 ! snow albedo: update snow age 1019 IF (snow(i) .GT. 0.0001) THEN 1020 agesno(i) = (agesno(i) + (1.-agesno(i)/50.)*dtime/86400.) & 1021 *EXP(-1.*MAX(0.0, precip_snow(i))*dtime/5.) 1022 ELSE 1023 agesno(i) = 0.0 1024 END IF 1025 ! snow albedo 1026 alb_snow = alb_sno_wet + (alb_sno_dry - alb_sno_wet)*EXP(-agesno(i)/50.) 1027 ! ice albedo (varies with ice thickness and temp) 1028 alb_ice = MAX(0.0, 0.13*LOG(100.*(hice(ki)+0.001)) + 0.1) 1029 !alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1) 1030 IF (tice(ki) .GT. t_freeze - 0.01) THEN 1031 alb_ice = MIN(alb_ice, alb_ice_wet) 1032 ELSE 1033 alb_ice = MIN(alb_ice, alb_ice_dry) 1034 END IF 1035 ! pond albedo 1036 alb_pond = 0.36 - 0.1*(2.0 + MIN(0.0, MAX(tice(ki) - t_melt, -2.0))) 1037 ! pond fraction 1038 frac_pond = 0.2*(2.0 + MIN(0.0, MAX(tice(ki) - t_melt, -2.0))) 1039 ! snow fraction 1040 frac_snow = MAX(0.0, MIN(1.0 - frac_pond, snow(i)/snow_min)) 1041 ! ice fraction 1042 frac_ice = MAX(0.0, 1.-frac_pond - frac_snow) 1043 ! total albedo 1044 alb1_new(i) = alb_snow*frac_snow + alb_ice*frac_ice + alb_pond*frac_pond 1045 END DO 1046 alb2_new(:) = alb1_new(:) 1047 1048 END SELECT 1049 1050 ! ------ End Albedo ---------- 1051 1052 END SUBROUTINE ocean_forced_ice 1053 !************************************************************************************* 1054 1055 SUBROUTINE ocean_forced_ice_reset_bilg_cumul(itime, dtime, bilg_cumul) 1056 USE dimphy, ONLY: klon 1057 USE surface_data, ONLY: iflag_seaice, type_ocean, version_ocean 1058 IMPLICIT NONE 1059 INTEGER, INTENT(IN) :: itime 1060 REAL, INTENT(IN) :: dtime 1061 REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul 1062 INTEGER :: cpl_pas 1063 1064 IF (iflag_seaice /= 0) THEN 1065 cpl_pas = NINT(86400./dtime*1.0) ! une fois par jour 1066 IF (MOD(itime, cpl_pas) .EQ. 0) THEN ! time to update tslab 1067 IF (.NOT. (type_ocean == 'couple' .OR. (type_ocean == 'slab' .AND. version_ocean == 'sicINT'))) THEN 1068 bilg_cumul(1:klon) = 0. 1069 END IF 1070 END IF ! coupling time 1071 END IF 1072 1073 END SUBROUTINE ocean_forced_ice_reset_bilg_cumul 1074 1133 1075 END MODULE ocean_forced_mod 1134 1076 1135 1136 1137 1138 1139 -
LMDZ6/trunk/libf/phylmd/phyetat0_mod.f90
r6059 r6070 18 18 !GG USE surface_data, ONLY : type_ocean, version_ocean 19 19 USE surface_data, ONLY : type_ocean, version_ocean, iflag_seaice, & 20 iflag_seaice_alb, iflag_leads 20 iflag_seaice_alb, iflag_leads, sic_hice_fixed 21 21 !GG 22 22 USE phyetat0_get_mod, ONLY : phyetat0_get, phyetat0_srf … … 725 725 !IF (iflag_seaice == 2) THEN 726 726 727 ! initial ice thickness (if not read in the startphy) could be included in the .def files 727 728 found=phyetat0_get(hice,"hice","Ice thickness",0.) 728 729 IF (.NOT. found) THEN 729 730 PRINT*, "phyetat0: Le champ <hice> est absent" 730 731 PRINT*, "Initialisation a hice=1m " 731 hice(:)= 1.0732 hice(:)=sic_hice_fixed 732 733 END IF 733 734 found=phyetat0_get(tice,"tice","Sea Ice temperature",0.) -
LMDZ6/trunk/libf/phylmd/surf_ocean_mod.F90
r5989 r6070 21 21 tsurf_new, dflux_s, dflux_l, lmt_bils, & 22 22 flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, & 23 !GG dt_ds, tkt, tks, taur, sss)24 23 dt_ds, tkt, tks, taur, sss, & 25 24 dthetadz300,Ampl & 26 !GG27 25 #ifdef ISO 28 26 & ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, & … … 282 280 radsol, snow, agesno, & 283 281 qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & 284 !GG tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa)285 282 tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa, & 286 283 dthetadz300, pctsrf, Ampl & 287 !GG288 284 #ifdef ISO 289 285 ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, & -
LMDZ6/trunk/libf/phylmd/surface_data.f90
r5971 r6070 59 59 REAL,SAVE :: fseaS 60 60 !$OMP THREADPRIVATE(fseaS) 61 !GG 61 REAL,SAVE :: sic_hice_fixed 62 !$OMP THREADPRIVATE(sic_hice_fixed) 63 62 64 63 65 ! if type_ocean=couple : version_ocean=opa8 ou nemo
Note: See TracChangeset
for help on using the changeset viewer.
