| 1 | !Completed |
|---|
| 2 | MODULE ocean_slab_mod |
|---|
| 3 | |
|---|
| 4 | ! This module is used for both surface ocean and sea-ice when using the slab ocean, |
|---|
| 5 | ! "ocean=slab". |
|---|
| 6 | |
|---|
| 7 | USE dimphy |
|---|
| 8 | USE indice_sol_mod |
|---|
| 9 | USE surface_data |
|---|
| 10 | USE lmdz_grid_phy, ONLY: klon_glo |
|---|
| 11 | USE lmdz_phys_mpi_data, ONLY: is_mpi_root |
|---|
| 12 | USE lmdz_abort_physic, ONLY: abort_physic |
|---|
| 13 | |
|---|
| 14 | IMPLICIT NONE |
|---|
| 15 | PRIVATE |
|---|
| 16 | PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice |
|---|
| 17 | |
|---|
| 18 | !*********************************************************************************** |
|---|
| 19 | ! Global saved variables |
|---|
| 20 | !*********************************************************************************** |
|---|
| 21 | ! number of slab vertical layers |
|---|
| 22 | INTEGER, PUBLIC, SAVE :: nslay |
|---|
| 23 | !$OMP THREADPRIVATE(nslay) |
|---|
| 24 | ! timestep for coupling (update slab temperature) in timesteps |
|---|
| 25 | INTEGER, PRIVATE, SAVE :: cpl_pas |
|---|
| 26 | !$OMP THREADPRIVATE(cpl_pas) |
|---|
| 27 | ! cyang = 1/heat capacity of top layer (rho.c.H) |
|---|
| 28 | REAL, PRIVATE, SAVE :: cyang |
|---|
| 29 | !$OMP THREADPRIVATE(cyang) |
|---|
| 30 | ! depth of slab layers (1 or 2) |
|---|
| 31 | REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh |
|---|
| 32 | !$OMP THREADPRIVATE(slabh) |
|---|
| 33 | ! slab temperature |
|---|
| 34 | REAL, ALLOCATABLE, DIMENSION(:, :), PUBLIC, SAVE :: tslab |
|---|
| 35 | !$OMP THREADPRIVATE(tslab) |
|---|
| 36 | ! heat flux convergence due to Ekman |
|---|
| 37 | REAL, ALLOCATABLE, DIMENSION(:, :), PUBLIC, SAVE :: dt_ekman |
|---|
| 38 | !$OMP THREADPRIVATE(dt_ekman) |
|---|
| 39 | ! heat flux convergence due to horiz diffusion |
|---|
| 40 | REAL, ALLOCATABLE, DIMENSION(:, :), PUBLIC, SAVE :: dt_hdiff |
|---|
| 41 | !$OMP THREADPRIVATE(dt_hdiff) |
|---|
| 42 | ! heat flux convergence due to GM eddy advection |
|---|
| 43 | REAL, ALLOCATABLE, DIMENSION(:, :), PUBLIC, SAVE :: dt_gm |
|---|
| 44 | !$OMP THREADPRIVATE(dt_gm) |
|---|
| 45 | ! Heat Flux correction |
|---|
| 46 | REAL, ALLOCATABLE, DIMENSION(:, :), PUBLIC, SAVE :: dt_qflux |
|---|
| 47 | !$OMP THREADPRIVATE(dt_qflux) |
|---|
| 48 | ! fraction of ocean covered by sea ice (sic / (oce+sic)) |
|---|
| 49 | REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic |
|---|
| 50 | !$OMP THREADPRIVATE(fsic) |
|---|
| 51 | ! temperature of the sea ice |
|---|
| 52 | REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice |
|---|
| 53 | !$OMP THREADPRIVATE(tice) |
|---|
| 54 | ! sea ice thickness, in kg/m2 |
|---|
| 55 | REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice |
|---|
| 56 | !$OMP THREADPRIVATE(seaice) |
|---|
| 57 | ! net surface heat flux, weighted by open ocean fraction |
|---|
| 58 | ! slab_bils accumulated over cpl_pas timesteps |
|---|
| 59 | REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum |
|---|
| 60 | !$OMP THREADPRIVATE(bils_cum) |
|---|
| 61 | ! net heat flux into the ocean below the ice : conduction + solar radiation |
|---|
| 62 | REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg |
|---|
| 63 | !$OMP THREADPRIVATE(slab_bilg) |
|---|
| 64 | ! slab_bilg over cpl_pas timesteps |
|---|
| 65 | REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum |
|---|
| 66 | !$OMP THREADPRIVATE(bilg_cum) |
|---|
| 67 | ! wind stress saved over cpl_pas timesteps |
|---|
| 68 | REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: taux_cum |
|---|
| 69 | !$OMP THREADPRIVATE(taux_cum) |
|---|
| 70 | REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tauy_cum |
|---|
| 71 | !$OMP THREADPRIVATE(tauy_cum) |
|---|
| 72 | |
|---|
| 73 | !*********************************************************************************** |
|---|
| 74 | ! Parameters (could be read in def file: move to slab_init) |
|---|
| 75 | !*********************************************************************************** |
|---|
| 76 | ! snow and ice physical characteristics: |
|---|
| 77 | REAL, PARAMETER :: t_freeze = 271.35 ! freezing sea water temp |
|---|
| 78 | REAL, PARAMETER :: t_melt = 273.15 ! melting ice temp |
|---|
| 79 | REAL, PARAMETER :: sno_den = 300. !mean snow density, kg/m3 |
|---|
| 80 | REAL, PARAMETER :: ice_den = 917. ! ice density |
|---|
| 81 | REAL, PARAMETER :: sea_den = 1025. ! sea water density |
|---|
| 82 | REAL, PARAMETER :: ice_cond = 2.17 * ice_den !conductivity of ice |
|---|
| 83 | REAL, PARAMETER :: sno_cond = 0.31 * sno_den ! conductivity of snow |
|---|
| 84 | REAL, PARAMETER :: ice_cap = 2067. ! specific heat capacity, snow and ice |
|---|
| 85 | REAL, PARAMETER :: sea_cap = 3995. ! specific heat capacity, snow and ice |
|---|
| 86 | REAL, PARAMETER :: ice_lat = 334000. ! freeze /melt latent heat snow and ice |
|---|
| 87 | |
|---|
| 88 | ! control of snow and ice cover & freeze / melt (heights converted to kg/m2) |
|---|
| 89 | REAL, PARAMETER :: snow_min = 0.05 * sno_den !critical snow height 5 cm |
|---|
| 90 | REAL, PARAMETER :: snow_wfact = 0.4 ! max fraction of falling snow blown into ocean |
|---|
| 91 | REAL, PARAMETER :: ice_frac_min = 0.001 |
|---|
| 92 | REAL, PARAMETER :: ice_frac_max = 1. ! less than 1. if min leads fraction |
|---|
| 93 | REAL, PARAMETER :: h_ice_min = 0.01 * ice_den ! min ice thickness |
|---|
| 94 | REAL, PARAMETER :: h_ice_thin = 0.15 * ice_den ! thin ice thickness |
|---|
| 95 | ! below ice_thin, priority is melt lateral / grow height |
|---|
| 96 | ! ice_thin is also height of new ice |
|---|
| 97 | REAL, PARAMETER :: h_ice_thick = 2.5 * ice_den ! thin ice thickness |
|---|
| 98 | ! above ice_thick, priority is melt height / grow lateral |
|---|
| 99 | REAL, PARAMETER :: h_ice_new = 1. * ice_den ! max height of new open ocean ice |
|---|
| 100 | REAL, PARAMETER :: h_ice_max = 10. * ice_den ! max ice height |
|---|
| 101 | |
|---|
| 102 | ! albedo and radiation parameters |
|---|
| 103 | REAL, PARAMETER :: alb_sno_min = 0.55 !min snow albedo |
|---|
| 104 | REAL, PARAMETER :: alb_sno_del = 0.3 !max snow albedo = min + del |
|---|
| 105 | REAL, PARAMETER :: alb_ice_dry = 0.75 !dry thick ice |
|---|
| 106 | REAL, PARAMETER :: alb_ice_wet = 0.66 !melting thick ice |
|---|
| 107 | REAL, PARAMETER :: pen_frac = 0.3 !fraction of shortwave penetrating into the |
|---|
| 108 | ! ice (no snow) |
|---|
| 109 | REAL, PARAMETER :: pen_ext = 1.5 !extinction of penetrating shortwave (m-1) |
|---|
| 110 | |
|---|
| 111 | ! horizontal transport |
|---|
| 112 | LOGICAL, PUBLIC, SAVE :: slab_hdiff |
|---|
| 113 | !$OMP THREADPRIVATE(slab_hdiff) |
|---|
| 114 | LOGICAL, PUBLIC, SAVE :: slab_gm |
|---|
| 115 | !$OMP THREADPRIVATE(slab_gm) |
|---|
| 116 | REAL, PRIVATE, SAVE :: coef_hdiff ! coefficient for horizontal diffusion |
|---|
| 117 | !$OMP THREADPRIVATE(coef_hdiff) |
|---|
| 118 | INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment |
|---|
| 119 | !$OMP THREADPRIVATE(slab_ekman) |
|---|
| 120 | !$OMP THREADPRIVATE(slab_cadj) |
|---|
| 121 | |
|---|
| 122 | !*********************************************************************************** |
|---|
| 123 | |
|---|
| 124 | CONTAINS |
|---|
| 125 | |
|---|
| 126 | !*********************************************************************************** |
|---|
| 127 | |
|---|
| 128 | SUBROUTINE ocean_slab_init(dtime, pctsrf_rst) |
|---|
| 129 | !, seaice_rst etc |
|---|
| 130 | |
|---|
| 131 | USE lmdz_ioipsl_getin_p, ONLY: getin_p |
|---|
| 132 | USE lmdz_phys_transfert_para, ONLY: gather |
|---|
| 133 | USE slab_heat_transp_mod, ONLY: ini_slab_transp |
|---|
| 134 | |
|---|
| 135 | ! Input variables |
|---|
| 136 | !*********************************************************************************** |
|---|
| 137 | REAL, INTENT(IN) :: dtime |
|---|
| 138 | ! Variables read from restart file |
|---|
| 139 | REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst |
|---|
| 140 | ! surface fractions from start file |
|---|
| 141 | |
|---|
| 142 | ! Local variables |
|---|
| 143 | !************************************************************************************ |
|---|
| 144 | INTEGER :: error |
|---|
| 145 | REAL, DIMENSION(klon_glo) :: zmasq_glo |
|---|
| 146 | CHARACTER (len = 80) :: abort_message |
|---|
| 147 | CHARACTER (len = 20) :: modname = 'ocean_slab_intit' |
|---|
| 148 | |
|---|
| 149 | !*********************************************************************************** |
|---|
| 150 | ! Define some parameters |
|---|
| 151 | !*********************************************************************************** |
|---|
| 152 | ! Number of slab layers |
|---|
| 153 | nslay = 2 |
|---|
| 154 | CALL getin_p('slab_layers', nslay) |
|---|
| 155 | print *, 'number of slab layers : ', nslay |
|---|
| 156 | ! Layer thickness |
|---|
| 157 | ALLOCATE(slabh(nslay), stat = error) |
|---|
| 158 | IF (error /= 0) THEN |
|---|
| 159 | abort_message = 'Pb allocation slabh' |
|---|
| 160 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 161 | ENDIF |
|---|
| 162 | slabh(1) = 50. |
|---|
| 163 | CALL getin_p('slab_depth', slabh(1)) |
|---|
| 164 | IF (nslay>1) THEN |
|---|
| 165 | slabh(2) = 150. |
|---|
| 166 | END IF |
|---|
| 167 | |
|---|
| 168 | ! cyang = 1/heat capacity of top layer (rho.c.H) |
|---|
| 169 | cyang = 1 / (slabh(1) * sea_den * sea_cap) |
|---|
| 170 | |
|---|
| 171 | ! cpl_pas coupling period (update of tslab and ice fraction) |
|---|
| 172 | ! pour un calcul a chaque pas de temps, cpl_pas=1 |
|---|
| 173 | cpl_pas = NINT(86400. / dtime * 1.0) ! une fois par jour |
|---|
| 174 | CALL getin_p('cpl_pas', cpl_pas) |
|---|
| 175 | print *, 'cpl_pas', cpl_pas |
|---|
| 176 | |
|---|
| 177 | ! Horizontal diffusion |
|---|
| 178 | slab_hdiff = .FALSE. |
|---|
| 179 | CALL getin_p('slab_hdiff', slab_hdiff) |
|---|
| 180 | coef_hdiff = 25000. |
|---|
| 181 | CALL getin_p('coef_hdiff', coef_hdiff) |
|---|
| 182 | ! Ekman transport |
|---|
| 183 | slab_ekman = 0 |
|---|
| 184 | CALL getin_p('slab_ekman', slab_ekman) |
|---|
| 185 | ! GM eddy advection (2-layers only) |
|---|
| 186 | slab_gm = .FALSE. |
|---|
| 187 | CALL getin_p('slab_gm', slab_gm) |
|---|
| 188 | IF (slab_ekman<2) THEN |
|---|
| 189 | slab_gm = .FALSE. |
|---|
| 190 | ENDIF |
|---|
| 191 | ! Convective adjustment |
|---|
| 192 | IF (nslay==1) THEN |
|---|
| 193 | slab_cadj = 0 |
|---|
| 194 | ELSE |
|---|
| 195 | slab_cadj = 1 |
|---|
| 196 | END IF |
|---|
| 197 | CALL getin_p('slab_cadj', slab_cadj) |
|---|
| 198 | |
|---|
| 199 | !************************************************************************************ |
|---|
| 200 | ! Allocate surface fraction read from restart file |
|---|
| 201 | !************************************************************************************ |
|---|
| 202 | ALLOCATE(fsic(klon), stat = error) |
|---|
| 203 | IF (error /= 0) THEN |
|---|
| 204 | abort_message = 'Pb allocation tmp_pctsrf_slab' |
|---|
| 205 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 206 | ENDIF |
|---|
| 207 | fsic(:) = 0. |
|---|
| 208 | !zmasq = continent fraction |
|---|
| 209 | WHERE (1. - zmasq(:)>EPSFRA) |
|---|
| 210 | fsic(:) = pctsrf_rst(:, is_sic) / (1. - zmasq(:)) |
|---|
| 211 | END WHERE |
|---|
| 212 | |
|---|
| 213 | !************************************************************************************ |
|---|
| 214 | ! Allocate saved fields |
|---|
| 215 | !************************************************************************************ |
|---|
| 216 | ALLOCATE(tslab(klon, nslay), stat = error) |
|---|
| 217 | IF (error /= 0) CALL abort_physic & |
|---|
| 218 | (modname, 'pb allocation tslab', 1) |
|---|
| 219 | |
|---|
| 220 | ALLOCATE(bils_cum(klon), stat = error) |
|---|
| 221 | IF (error /= 0) THEN |
|---|
| 222 | abort_message = 'Pb allocation slab_bils_cum' |
|---|
| 223 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 224 | ENDIF |
|---|
| 225 | bils_cum(:) = 0.0 |
|---|
| 226 | |
|---|
| 227 | IF (version_ocean=='sicINT') THEN ! interactive sea ice |
|---|
| 228 | ALLOCATE(slab_bilg(klon), stat = error) |
|---|
| 229 | IF (error /= 0) THEN |
|---|
| 230 | abort_message = 'Pb allocation slab_bilg' |
|---|
| 231 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 232 | ENDIF |
|---|
| 233 | slab_bilg(:) = 0.0 |
|---|
| 234 | ALLOCATE(bilg_cum(klon), stat = error) |
|---|
| 235 | IF (error /= 0) THEN |
|---|
| 236 | abort_message = 'Pb allocation slab_bilg_cum' |
|---|
| 237 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 238 | ENDIF |
|---|
| 239 | bilg_cum(:) = 0.0 |
|---|
| 240 | ALLOCATE(tice(klon), stat = error) |
|---|
| 241 | IF (error /= 0) THEN |
|---|
| 242 | abort_message = 'Pb allocation slab_tice' |
|---|
| 243 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 244 | ENDIF |
|---|
| 245 | ALLOCATE(seaice(klon), stat = error) |
|---|
| 246 | IF (error /= 0) THEN |
|---|
| 247 | abort_message = 'Pb allocation slab_seaice' |
|---|
| 248 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 249 | ENDIF |
|---|
| 250 | END IF |
|---|
| 251 | |
|---|
| 252 | IF (slab_hdiff) THEN !horizontal diffusion |
|---|
| 253 | ALLOCATE(dt_hdiff(klon, nslay), stat = error) |
|---|
| 254 | IF (error /= 0) THEN |
|---|
| 255 | abort_message = 'Pb allocation dt_hdiff' |
|---|
| 256 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 257 | ENDIF |
|---|
| 258 | dt_hdiff(:, :) = 0.0 |
|---|
| 259 | ENDIF |
|---|
| 260 | |
|---|
| 261 | ALLOCATE(dt_qflux(klon, nslay), stat = error) |
|---|
| 262 | IF (error /= 0) THEN |
|---|
| 263 | abort_message = 'Pb allocation dt_qflux' |
|---|
| 264 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 265 | ENDIF |
|---|
| 266 | dt_qflux(:, :) = 0.0 |
|---|
| 267 | |
|---|
| 268 | IF (slab_gm) THEN !GM advection |
|---|
| 269 | ALLOCATE(dt_gm(klon, nslay), stat = error) |
|---|
| 270 | IF (error /= 0) THEN |
|---|
| 271 | abort_message = 'Pb allocation dt_gm' |
|---|
| 272 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 273 | ENDIF |
|---|
| 274 | dt_gm(:, :) = 0.0 |
|---|
| 275 | ENDIF |
|---|
| 276 | |
|---|
| 277 | IF (slab_ekman>0) THEN ! ekman transport |
|---|
| 278 | ALLOCATE(dt_ekman(klon, nslay), stat = error) |
|---|
| 279 | IF (error /= 0) THEN |
|---|
| 280 | abort_message = 'Pb allocation dt_ekman' |
|---|
| 281 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 282 | ENDIF |
|---|
| 283 | dt_ekman(:, :) = 0.0 |
|---|
| 284 | ALLOCATE(taux_cum(klon), stat = error) |
|---|
| 285 | IF (error /= 0) THEN |
|---|
| 286 | abort_message = 'Pb allocation taux_cum' |
|---|
| 287 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 288 | ENDIF |
|---|
| 289 | taux_cum(:) = 0.0 |
|---|
| 290 | ALLOCATE(tauy_cum(klon), stat = error) |
|---|
| 291 | IF (error /= 0) THEN |
|---|
| 292 | abort_message = 'Pb allocation tauy_cum' |
|---|
| 293 | CALL abort_physic(modname, abort_message, 1) |
|---|
| 294 | ENDIF |
|---|
| 295 | tauy_cum(:) = 0.0 |
|---|
| 296 | ENDIF |
|---|
| 297 | |
|---|
| 298 | ! Initialize transport |
|---|
| 299 | IF (slab_hdiff.OR.(slab_ekman>0)) THEN |
|---|
| 300 | CALL gather(zmasq, zmasq_glo) |
|---|
| 301 | ! Master thread/process only |
|---|
| 302 | !$OMP MASTER |
|---|
| 303 | IF (is_mpi_root) THEN |
|---|
| 304 | CALL ini_slab_transp(zmasq_glo) |
|---|
| 305 | END IF |
|---|
| 306 | !$OMP END MASTER |
|---|
| 307 | END IF |
|---|
| 308 | |
|---|
| 309 | END SUBROUTINE ocean_slab_init |
|---|
| 310 | |
|---|
| 311 | !*********************************************************************************** |
|---|
| 312 | |
|---|
| 313 | SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified) |
|---|
| 314 | |
|---|
| 315 | ! this routine sends back the sea ice and ocean fraction to the main physics |
|---|
| 316 | ! routine. Called only with interactive sea ice |
|---|
| 317 | |
|---|
| 318 | ! Arguments |
|---|
| 319 | !************************************************************************************ |
|---|
| 320 | INTEGER, INTENT(IN) :: itime ! current timestep |
|---|
| 321 | INTEGER, INTENT(IN) :: jour ! day in year (not |
|---|
| 322 | REAL, INTENT(IN) :: dtime ! physics timestep (s) |
|---|
| 323 | REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: pctsrf_chg ! sub-surface fraction |
|---|
| 324 | LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is |
|---|
| 325 | ! modified at this time step |
|---|
| 326 | |
|---|
| 327 | pctsrf_chg(:, is_oce) = (1. - fsic(:)) * (1. - zmasq(:)) |
|---|
| 328 | pctsrf_chg(:, is_sic) = fsic(:) * (1. - zmasq(:)) |
|---|
| 329 | is_modified = .TRUE. |
|---|
| 330 | |
|---|
| 331 | END SUBROUTINE ocean_slab_frac |
|---|
| 332 | |
|---|
| 333 | !************************************************************************************ |
|---|
| 334 | |
|---|
| 335 | SUBROUTINE ocean_slab_noice(& |
|---|
| 336 | itime, dtime, jour, knon, knindex, & |
|---|
| 337 | p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, & |
|---|
| 338 | AcoefH, AcoefQ, BcoefH, BcoefQ, & |
|---|
| 339 | AcoefU, AcoefV, BcoefU, BcoefV, & |
|---|
| 340 | ps, u1, v1, gustiness, tsurf_in, & |
|---|
| 341 | radsol, snow, & |
|---|
| 342 | qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & |
|---|
| 343 | tsurf_new, dflux_s, dflux_l, slab_bils) |
|---|
| 344 | |
|---|
| 345 | USE calcul_fluxs_mod |
|---|
| 346 | USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2, slab_gmdiff |
|---|
| 347 | USE lmdz_phys_para |
|---|
| 348 | |
|---|
| 349 | INCLUDE "clesphys.h" |
|---|
| 350 | |
|---|
| 351 | ! This routine |
|---|
| 352 | ! (1) computes surface turbulent fluxes over points with some open ocean |
|---|
| 353 | ! (2) reads additional Q-flux (everywhere) |
|---|
| 354 | ! (3) computes horizontal transport (diffusion & Ekman) |
|---|
| 355 | ! (4) updates slab temperature every cpl_pas ; creates new ice if needed. |
|---|
| 356 | |
|---|
| 357 | ! Note : |
|---|
| 358 | ! klon total number of points |
|---|
| 359 | ! knon number of points with open ocean (varies with time) |
|---|
| 360 | ! knindex gives position of the knon points within klon. |
|---|
| 361 | ! In general, local saved variables have klon values |
|---|
| 362 | ! variables exchanged with PBL module have knon. |
|---|
| 363 | |
|---|
| 364 | ! Input arguments |
|---|
| 365 | !*********************************************************************************** |
|---|
| 366 | INTEGER, INTENT(IN) :: itime ! current timestep INTEGER, |
|---|
| 367 | INTEGER, INTENT(IN) :: jour ! day in year (for Q-Flux) |
|---|
| 368 | INTEGER, INTENT(IN) :: knon ! number of points |
|---|
| 369 | INTEGER, DIMENSION(klon), INTENT(IN) :: knindex |
|---|
| 370 | REAL, INTENT(IN) :: dtime ! timestep (s) |
|---|
| 371 | REAL, DIMENSION(klon), INTENT(IN) :: p1lay |
|---|
| 372 | REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm |
|---|
| 373 | ! drag coefficients |
|---|
| 374 | REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow |
|---|
| 375 | REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum ! near surface T, q |
|---|
| 376 | REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ |
|---|
| 377 | REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV |
|---|
| 378 | ! exchange coefficients for boundary layer scheme |
|---|
| 379 | REAL, DIMENSION(klon), INTENT(IN) :: ps ! surface pressure |
|---|
| 380 | REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness ! surface wind |
|---|
| 381 | REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in ! surface temperature |
|---|
| 382 | REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux |
|---|
| 383 | |
|---|
| 384 | ! In/Output arguments |
|---|
| 385 | !************************************************************************************ |
|---|
| 386 | REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2 |
|---|
| 387 | |
|---|
| 388 | ! Output arguments |
|---|
| 389 | !************************************************************************************ |
|---|
| 390 | REAL, DIMENSION(klon), INTENT(OUT) :: qsurf |
|---|
| 391 | REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat |
|---|
| 392 | REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 |
|---|
| 393 | REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! new surface tempearture |
|---|
| 394 | REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l |
|---|
| 395 | REAL, DIMENSION(klon), INTENT(OUT) :: slab_bils |
|---|
| 396 | |
|---|
| 397 | ! Local variables |
|---|
| 398 | !************************************************************************************ |
|---|
| 399 | INTEGER :: i, ki, k |
|---|
| 400 | REAL :: t_cadj |
|---|
| 401 | ! for surface heat fluxes |
|---|
| 402 | REAL, DIMENSION(klon) :: cal, beta, dif_grnd |
|---|
| 403 | ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes |
|---|
| 404 | REAL, DIMENSION(klon) :: diff_sst, diff_siv |
|---|
| 405 | REAL, DIMENSION(klon, nslay) :: lmt_bils |
|---|
| 406 | ! for surface wind stress |
|---|
| 407 | REAL, DIMENSION(klon) :: u0, v0 |
|---|
| 408 | REAL, DIMENSION(klon) :: u1_lay, v1_lay |
|---|
| 409 | ! for new ice creation |
|---|
| 410 | REAL :: e_freeze, h_new, dfsic |
|---|
| 411 | ! horizontal diffusion and Ekman local vars |
|---|
| 412 | ! dimension = global domain (klon_glo) instead of // subdomains |
|---|
| 413 | REAL, DIMENSION(klon_glo, nslay) :: dt_hdiff_glo, dt_ekman_glo, dt_gm_glo |
|---|
| 414 | ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop |
|---|
| 415 | REAL, DIMENSION(klon_glo, nslay) :: dt_hdiff_tmp, dt_ekman_tmp |
|---|
| 416 | REAL, DIMENSION(klon_glo, nslay) :: tslab_glo |
|---|
| 417 | REAL, DIMENSION(klon_glo) :: taux_glo, tauy_glo |
|---|
| 418 | |
|---|
| 419 | !**************************************************************************************** |
|---|
| 420 | ! 1) Surface fluxes calculation |
|---|
| 421 | |
|---|
| 422 | !**************************************************************************************** |
|---|
| 423 | !cal(:) = 0. ! infinite thermal inertia |
|---|
| 424 | !beta(:) = 1. ! wet surface |
|---|
| 425 | !dif_grnd(:) = 0. ! no diffusion into ground |
|---|
| 426 | ! EV: use calbeta |
|---|
| 427 | CALL calbeta(dtime, is_oce, knon, snow, qsurf, beta, cal, dif_grnd) |
|---|
| 428 | |
|---|
| 429 | |
|---|
| 430 | |
|---|
| 431 | ! Suppose zero surface speed |
|---|
| 432 | u0(:) = 0.0 |
|---|
| 433 | v0(:) = 0.0 |
|---|
| 434 | u1_lay(:) = u1(:) - u0(:) |
|---|
| 435 | v1_lay(:) = v1(:) - v0(:) |
|---|
| 436 | |
|---|
| 437 | ! Compute latent & sensible fluxes |
|---|
| 438 | CALL calcul_fluxs(knon, is_oce, dtime, & |
|---|
| 439 | tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, & |
|---|
| 440 | precip_rain, precip_snow, snow, qsurf, & |
|---|
| 441 | radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & |
|---|
| 442 | f_qsat_oce, AcoefH, AcoefQ, BcoefH, BcoefQ, & |
|---|
| 443 | tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) |
|---|
| 444 | |
|---|
| 445 | ! save total cumulated heat fluxes locally |
|---|
| 446 | ! radiative + turbulent + melt of falling snow |
|---|
| 447 | slab_bils(:) = 0. |
|---|
| 448 | DO i = 1, knon |
|---|
| 449 | ki = knindex(i) |
|---|
| 450 | slab_bils(ki) = (1. - fsic(ki)) * (fluxlat(i) + fluxsens(i) + radsol(i) & |
|---|
| 451 | - precip_snow(i) * ice_lat * (1. + snow_wfact * fsic(ki))) |
|---|
| 452 | bils_cum(ki) = bils_cum(ki) + slab_bils(ki) |
|---|
| 453 | END DO |
|---|
| 454 | |
|---|
| 455 | ! Compute surface wind stress |
|---|
| 456 | CALL calcul_flux_wind(knon, dtime, & |
|---|
| 457 | u0, v0, u1, v1, gustiness, cdragm, & |
|---|
| 458 | AcoefU, AcoefV, BcoefU, BcoefV, & |
|---|
| 459 | p1lay, temp_air, & |
|---|
| 460 | flux_u1, flux_v1) |
|---|
| 461 | |
|---|
| 462 | ! save cumulated wind stress |
|---|
| 463 | IF (slab_ekman>0) THEN |
|---|
| 464 | DO i = 1, knon |
|---|
| 465 | ki = knindex(i) |
|---|
| 466 | taux_cum(ki) = taux_cum(ki) + flux_u1(i) * (1. - fsic(ki)) / cpl_pas |
|---|
| 467 | tauy_cum(ki) = tauy_cum(ki) + flux_v1(i) * (1. - fsic(ki)) / cpl_pas |
|---|
| 468 | END DO |
|---|
| 469 | ENDIF |
|---|
| 470 | |
|---|
| 471 | !**************************************************************************************** |
|---|
| 472 | ! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc |
|---|
| 473 | |
|---|
| 474 | !**************************************************************************************** |
|---|
| 475 | CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) |
|---|
| 476 | ! lmt_bils and diff_sst,siv saved by limit_slab |
|---|
| 477 | ! qflux = total QFlux correction (in W/m2) |
|---|
| 478 | dt_qflux(:, 1) = lmt_bils(:, 1) + diff_sst(:) / cyang / 86400. - diff_siv(:) * ice_den * ice_lat / 86400. |
|---|
| 479 | IF (nslay>1) THEN |
|---|
| 480 | dt_qflux(:, 2:nslay) = lmt_bils(:, 2:nslay) |
|---|
| 481 | END IF |
|---|
| 482 | |
|---|
| 483 | !**************************************************************************************** |
|---|
| 484 | ! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport) |
|---|
| 485 | ! Bring to freezing temp and make sea ice if necessary |
|---|
| 486 | |
|---|
| 487 | !***********************************************o***************************************** |
|---|
| 488 | tsurf_new = tsurf_in |
|---|
| 489 | IF (MOD(itime, cpl_pas)==0) THEN ! time to update tslab & fraction |
|---|
| 490 | ! *********************************** |
|---|
| 491 | ! Horizontal transport |
|---|
| 492 | ! *********************************** |
|---|
| 493 | IF (slab_ekman>0) THEN |
|---|
| 494 | ! copy wind stress to global var |
|---|
| 495 | CALL gather(taux_cum, taux_glo) |
|---|
| 496 | CALL gather(tauy_cum, tauy_glo) |
|---|
| 497 | END IF |
|---|
| 498 | |
|---|
| 499 | IF (slab_hdiff.OR.(slab_ekman>0)) THEN |
|---|
| 500 | CALL gather(tslab, tslab_glo) |
|---|
| 501 | ! Compute horiz transport on one process only |
|---|
| 502 | IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus |
|---|
| 503 | IF (slab_hdiff) THEN |
|---|
| 504 | dt_hdiff_glo(:, :) = 0. |
|---|
| 505 | END IF |
|---|
| 506 | IF (slab_ekman>0) THEN |
|---|
| 507 | dt_ekman_glo(:, :) = 0. |
|---|
| 508 | END IF |
|---|
| 509 | IF (slab_gm) THEN |
|---|
| 510 | dt_gm_glo(:, :) = 0. |
|---|
| 511 | END IF |
|---|
| 512 | DO i = 1, cpl_pas ! time splitting for numerical stability |
|---|
| 513 | IF (slab_ekman>0) THEN |
|---|
| 514 | SELECT CASE (slab_ekman) |
|---|
| 515 | CASE (1) |
|---|
| 516 | CALL slab_ekman1(taux_glo, tauy_glo, tslab_glo, dt_ekman_tmp) |
|---|
| 517 | CASE (2) |
|---|
| 518 | CALL slab_ekman2(taux_glo, tauy_glo, tslab_glo, dt_ekman_tmp, dt_hdiff_tmp, slab_gm) |
|---|
| 519 | CASE DEFAULT |
|---|
| 520 | dt_ekman_tmp(:, :) = 0. |
|---|
| 521 | END SELECT |
|---|
| 522 | dt_ekman_glo(:, :) = dt_ekman_glo(:, :) + dt_ekman_tmp(:, :) |
|---|
| 523 | ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1 |
|---|
| 524 | DO k = 1, nslay |
|---|
| 525 | dt_ekman_tmp(:, k) = dt_ekman_tmp(:, k) / (slabh(k) * sea_den) |
|---|
| 526 | ENDDO |
|---|
| 527 | tslab_glo = tslab_glo + dt_ekman_tmp * dtime |
|---|
| 528 | IF (slab_gm) THEN ! Gent-McWilliams eddy advection |
|---|
| 529 | dt_gm_glo(:, :) = dt_gm_glo(:, :) + dt_hdiff_tmp(:, :) |
|---|
| 530 | ! convert dt from K.s-1.(kg.m-2) to K.s-1 |
|---|
| 531 | DO k = 1, nslay |
|---|
| 532 | dt_hdiff_tmp(:, k) = dt_hdiff_tmp(:, k) / (slabh(k) * sea_den) |
|---|
| 533 | END DO |
|---|
| 534 | tslab_glo = tslab_glo + dt_hdiff_tmp * dtime |
|---|
| 535 | END IF |
|---|
| 536 | ENDIF |
|---|
| 537 | ! GM included in Ekman_2 |
|---|
| 538 | ! IF (slab_gm) THEN ! Gent-McWilliams eddy advection |
|---|
| 539 | ! CALL slab_gmdiff(tslab_glo,dt_hdiff_tmp) |
|---|
| 540 | ! ! convert dt_gm from K.m.s-1 to K.s-1 |
|---|
| 541 | ! DO k=1,nslay |
|---|
| 542 | ! dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/slabh(k) |
|---|
| 543 | ! END DO |
|---|
| 544 | ! tslab_glo=tslab_glo+dt_hdiff_tmp*dtime |
|---|
| 545 | ! dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:) |
|---|
| 546 | ! END IF |
|---|
| 547 | IF (slab_hdiff) THEN ! horizontal diffusion |
|---|
| 548 | ! laplacian of slab T |
|---|
| 549 | CALL divgrad_phy(nslay, tslab_glo, dt_hdiff_tmp) |
|---|
| 550 | ! multiply by diff coef and normalize to 50m slab equivalent |
|---|
| 551 | dt_hdiff_tmp = dt_hdiff_tmp * coef_hdiff * 50. / SUM(slabh) |
|---|
| 552 | dt_hdiff_glo(:, :) = dt_hdiff_glo(:, :) + dt_hdiff_tmp(:, :) |
|---|
| 553 | tslab_glo = tslab_glo + dt_hdiff_tmp * dtime |
|---|
| 554 | END IF |
|---|
| 555 | END DO ! time splitting |
|---|
| 556 | IF (slab_hdiff) THEN |
|---|
| 557 | !dt_hdiff_glo saved in W/m2 |
|---|
| 558 | DO k = 1, nslay |
|---|
| 559 | dt_hdiff_glo(:, k) = dt_hdiff_glo(:, k) * slabh(k) * sea_den * sea_cap / cpl_pas |
|---|
| 560 | END DO |
|---|
| 561 | END IF |
|---|
| 562 | IF (slab_gm) THEN |
|---|
| 563 | !dt_hdiff_glo saved in W/m2 |
|---|
| 564 | dt_gm_glo(:, :) = dt_gm_glo(:, :) * sea_cap / cpl_pas |
|---|
| 565 | END IF |
|---|
| 566 | IF (slab_ekman>0) THEN |
|---|
| 567 | ! dt_ekman_glo saved in W/m2 |
|---|
| 568 | dt_ekman_glo(:, :) = dt_ekman_glo(:, :) * sea_cap / cpl_pas |
|---|
| 569 | END IF |
|---|
| 570 | END IF ! master process |
|---|
| 571 | !$OMP BARRIER |
|---|
| 572 | ! Send new fields back to all processes |
|---|
| 573 | CALL Scatter(tslab_glo, tslab) |
|---|
| 574 | IF (slab_hdiff) THEN |
|---|
| 575 | CALL Scatter(dt_hdiff_glo, dt_hdiff) |
|---|
| 576 | END IF |
|---|
| 577 | IF (slab_gm) THEN |
|---|
| 578 | CALL Scatter(dt_gm_glo, dt_gm) |
|---|
| 579 | END IF |
|---|
| 580 | IF (slab_ekman>0) THEN |
|---|
| 581 | CALL Scatter(dt_ekman_glo, dt_ekman) |
|---|
| 582 | ! clear wind stress |
|---|
| 583 | taux_cum(:) = 0. |
|---|
| 584 | tauy_cum(:) = 0. |
|---|
| 585 | END IF |
|---|
| 586 | ENDIF ! transport |
|---|
| 587 | |
|---|
| 588 | ! *********************************** |
|---|
| 589 | ! Other heat fluxes |
|---|
| 590 | ! *********************************** |
|---|
| 591 | ! Add read QFlux |
|---|
| 592 | DO k = 1, nslay |
|---|
| 593 | tslab(:, k) = tslab(:, k) + dt_qflux(:, k) * cyang * dtime * cpl_pas & |
|---|
| 594 | * slabh(1) / slabh(k) |
|---|
| 595 | END DO |
|---|
| 596 | ! Add cumulated surface fluxes |
|---|
| 597 | tslab(:, 1) = tslab(:, 1) + bils_cum(:) * cyang * dtime |
|---|
| 598 | ! Convective adjustment if 2 layers |
|---|
| 599 | IF ((nslay>1).AND.(slab_cadj>0)) THEN |
|---|
| 600 | DO i = 1, klon |
|---|
| 601 | IF (tslab(i, 2)>tslab(i, 1)) THEN |
|---|
| 602 | ! mean (mass-weighted) temperature |
|---|
| 603 | t_cadj = SUM(tslab(i, :) * slabh(:)) / SUM(slabh(:)) |
|---|
| 604 | tslab(i, 1) = t_cadj |
|---|
| 605 | tslab(i, 2) = t_cadj |
|---|
| 606 | END IF |
|---|
| 607 | END DO |
|---|
| 608 | END IF |
|---|
| 609 | ! *********************************** |
|---|
| 610 | ! Update surface temperature and ice |
|---|
| 611 | ! *********************************** |
|---|
| 612 | SELECT CASE(version_ocean) |
|---|
| 613 | CASE('sicNO') ! no sea ice even below freezing ! |
|---|
| 614 | DO i = 1, knon |
|---|
| 615 | ki = knindex(i) |
|---|
| 616 | tsurf_new(i) = tslab(ki, 1) |
|---|
| 617 | END DO |
|---|
| 618 | CASE('sicOBS') ! "realistic" case, for prescribed sea ice |
|---|
| 619 | ! tslab cannot be below freezing, or above it if there is sea ice |
|---|
| 620 | DO i = 1, knon |
|---|
| 621 | ki = knindex(i) |
|---|
| 622 | IF ((tslab(ki, 1)<t_freeze).OR.(fsic(ki)>epsfra)) THEN |
|---|
| 623 | tslab(ki, 1) = t_freeze |
|---|
| 624 | END IF |
|---|
| 625 | tsurf_new(i) = tslab(ki, 1) |
|---|
| 626 | END DO |
|---|
| 627 | CASE('sicINT') ! interactive sea ice |
|---|
| 628 | DO i = 1, knon |
|---|
| 629 | ki = knindex(i) |
|---|
| 630 | IF (fsic(ki)<epsfra) THEN ! Free of ice |
|---|
| 631 | IF (tslab(ki, 1)<t_freeze) THEN ! create new ice |
|---|
| 632 | ! quantity of new ice formed |
|---|
| 633 | e_freeze = (t_freeze - tslab(ki, 1)) / cyang / ice_lat |
|---|
| 634 | ! new ice |
|---|
| 635 | tice(ki) = t_freeze |
|---|
| 636 | fsic(ki) = MIN(ice_frac_max, e_freeze / h_ice_thin) |
|---|
| 637 | IF (fsic(ki)>ice_frac_min) THEN |
|---|
| 638 | seaice(ki) = MIN(e_freeze / fsic(ki), h_ice_max) |
|---|
| 639 | tslab(ki, 1) = t_freeze |
|---|
| 640 | ELSE |
|---|
| 641 | fsic(ki) = 0. |
|---|
| 642 | END IF |
|---|
| 643 | tsurf_new(i) = t_freeze |
|---|
| 644 | ELSE |
|---|
| 645 | tsurf_new(i) = tslab(ki, 1) |
|---|
| 646 | END IF |
|---|
| 647 | ELSE ! ice present |
|---|
| 648 | tsurf_new(i) = t_freeze |
|---|
| 649 | IF (tslab(ki, 1)<t_freeze) THEN ! create new ice |
|---|
| 650 | ! quantity of new ice formed over open ocean |
|---|
| 651 | e_freeze = (t_freeze - tslab(ki, 1)) / cyang * (1. - fsic(ki)) & |
|---|
| 652 | / (ice_lat + ice_cap / 2. * (t_freeze - tice(ki))) |
|---|
| 653 | ! new ice height and fraction |
|---|
| 654 | h_new = MIN(h_ice_new, seaice(ki)) ! max new height ice_new |
|---|
| 655 | dfsic = MIN(ice_frac_max - fsic(ki), e_freeze / h_new) |
|---|
| 656 | h_new = MIN(e_freeze / dfsic, h_ice_max) |
|---|
| 657 | ! update tslab to freezing over open ocean only |
|---|
| 658 | tslab(ki, 1) = tslab(ki, 1) * fsic(ki) + t_freeze * (1. - fsic(ki)) |
|---|
| 659 | ! update sea ice |
|---|
| 660 | seaice(ki) = (h_new * dfsic + seaice(ki) * fsic(ki)) & |
|---|
| 661 | / (dfsic + fsic(ki)) |
|---|
| 662 | fsic(ki) = fsic(ki) + dfsic |
|---|
| 663 | ! update snow? |
|---|
| 664 | END IF ! tslab below freezing |
|---|
| 665 | END IF ! sea ice present |
|---|
| 666 | END DO |
|---|
| 667 | END SELECT |
|---|
| 668 | bils_cum(:) = 0.0! clear cumulated fluxes |
|---|
| 669 | END IF ! coupling time |
|---|
| 670 | END SUBROUTINE ocean_slab_noice |
|---|
| 671 | |
|---|
| 672 | !**************************************************************************************** |
|---|
| 673 | |
|---|
| 674 | SUBROUTINE ocean_slab_ice(& |
|---|
| 675 | itime, dtime, jour, knon, knindex, & |
|---|
| 676 | tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, & |
|---|
| 677 | AcoefH, AcoefQ, BcoefH, BcoefQ, & |
|---|
| 678 | AcoefU, AcoefV, BcoefU, BcoefV, & |
|---|
| 679 | ps, u1, v1, gustiness, & |
|---|
| 680 | radsol, snow, qsurf, qsol, agesno, & |
|---|
| 681 | alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & |
|---|
| 682 | tsurf_new, dflux_s, dflux_l, swnet) |
|---|
| 683 | |
|---|
| 684 | USE calcul_fluxs_mod |
|---|
| 685 | |
|---|
| 686 | INCLUDE "YOMCST.h" |
|---|
| 687 | INCLUDE "clesphys.h" |
|---|
| 688 | |
|---|
| 689 | ! Input arguments |
|---|
| 690 | !**************************************************************************************** |
|---|
| 691 | INTEGER, INTENT(IN) :: itime, jour, knon |
|---|
| 692 | INTEGER, DIMENSION(klon), INTENT(IN) :: knindex |
|---|
| 693 | REAL, INTENT(IN) :: dtime |
|---|
| 694 | REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in |
|---|
| 695 | REAL, DIMENSION(klon), INTENT(IN) :: p1lay |
|---|
| 696 | REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm |
|---|
| 697 | REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow |
|---|
| 698 | REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum |
|---|
| 699 | REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ |
|---|
| 700 | REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV |
|---|
| 701 | REAL, DIMENSION(klon), INTENT(IN) :: ps |
|---|
| 702 | REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness |
|---|
| 703 | REAL, DIMENSION(klon), INTENT(IN) :: swnet |
|---|
| 704 | |
|---|
| 705 | ! In/Output arguments |
|---|
| 706 | !**************************************************************************************** |
|---|
| 707 | REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol |
|---|
| 708 | REAL, DIMENSION(klon), INTENT(INOUT) :: agesno |
|---|
| 709 | REAL, DIMENSION(klon), INTENT(INOUT) :: radsol |
|---|
| 710 | |
|---|
| 711 | ! Output arguments |
|---|
| 712 | !**************************************************************************************** |
|---|
| 713 | REAL, DIMENSION(klon), INTENT(OUT) :: qsurf |
|---|
| 714 | REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval |
|---|
| 715 | REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval |
|---|
| 716 | REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat |
|---|
| 717 | REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1 |
|---|
| 718 | REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new |
|---|
| 719 | REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l |
|---|
| 720 | |
|---|
| 721 | ! Local variables |
|---|
| 722 | !**************************************************************************************** |
|---|
| 723 | INTEGER :: i, ki |
|---|
| 724 | REAL, DIMENSION(klon) :: cal, beta, dif_grnd |
|---|
| 725 | REAL, DIMENSION(klon) :: u0, v0 |
|---|
| 726 | REAL, DIMENSION(klon) :: u1_lay, v1_lay |
|---|
| 727 | ! intermediate heat fluxes: |
|---|
| 728 | REAL :: f_cond, f_swpen |
|---|
| 729 | ! for snow/ice albedo: |
|---|
| 730 | REAL :: alb_snow, alb_ice, alb_pond |
|---|
| 731 | REAL :: frac_snow, frac_ice, frac_pond |
|---|
| 732 | ! for ice melt / freeze |
|---|
| 733 | REAL :: e_melt, snow_evap, h_test |
|---|
| 734 | ! dhsic, dfsic change in ice mass, fraction. |
|---|
| 735 | REAL :: dhsic, dfsic, frac_mf |
|---|
| 736 | |
|---|
| 737 | !**************************************************************************************** |
|---|
| 738 | ! 1) Flux calculation |
|---|
| 739 | !**************************************************************************************** |
|---|
| 740 | ! Suppose zero surface speed |
|---|
| 741 | u0(:) = 0.0 |
|---|
| 742 | v0(:) = 0.0 |
|---|
| 743 | u1_lay(:) = u1(:) - u0(:) |
|---|
| 744 | v1_lay(:) = v1(:) - v0(:) |
|---|
| 745 | |
|---|
| 746 | ! set beta, cal, compute conduction fluxes inside ice/snow |
|---|
| 747 | slab_bilg(:) = 0. |
|---|
| 748 | !dif_grnd(:)=0. |
|---|
| 749 | !beta(:) = 1. |
|---|
| 750 | ! EV: use calbeta to calculate beta and then recalculate properly cal |
|---|
| 751 | CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, cal, dif_grnd) |
|---|
| 752 | |
|---|
| 753 | DO i = 1, knon |
|---|
| 754 | ki = knindex(i) |
|---|
| 755 | IF (snow(i)>snow_min) THEN |
|---|
| 756 | ! snow-layer heat capacity |
|---|
| 757 | cal(i) = 2. * RCPD / (snow(i) * ice_cap) |
|---|
| 758 | ! snow conductive flux |
|---|
| 759 | f_cond = sno_cond * (tice(ki) - tsurf_in(i)) / snow(i) |
|---|
| 760 | ! all shortwave flux absorbed |
|---|
| 761 | f_swpen = 0. |
|---|
| 762 | ! bottom flux (ice conduction) |
|---|
| 763 | slab_bilg(ki) = ice_cond * (tice(ki) - t_freeze) / seaice(ki) |
|---|
| 764 | ! update ice temperature |
|---|
| 765 | tice(ki) = tice(ki) - 2. / ice_cap / (snow(i) + seaice(ki)) & |
|---|
| 766 | * (slab_bilg(ki) + f_cond) * dtime |
|---|
| 767 | ELSE ! bare ice |
|---|
| 768 | ! ice-layer heat capacity |
|---|
| 769 | cal(i) = 2. * RCPD / (seaice(ki) * ice_cap) |
|---|
| 770 | ! conductive flux |
|---|
| 771 | f_cond = ice_cond * (t_freeze - tice(ki)) / seaice(ki) |
|---|
| 772 | ! penetrative shortwave flux... |
|---|
| 773 | f_swpen = swnet(i) * pen_frac * exp(-pen_ext * seaice(ki) / ice_den) |
|---|
| 774 | slab_bilg(ki) = f_swpen - f_cond |
|---|
| 775 | END IF |
|---|
| 776 | radsol(i) = radsol(i) + f_cond - f_swpen |
|---|
| 777 | END DO |
|---|
| 778 | ! weight fluxes to ocean by sea ice fraction |
|---|
| 779 | slab_bilg(:) = slab_bilg(:) * fsic(:) |
|---|
| 780 | |
|---|
| 781 | ! calcul_fluxs (sens, lat etc) |
|---|
| 782 | CALL calcul_fluxs(knon, is_sic, dtime, & |
|---|
| 783 | tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, & |
|---|
| 784 | precip_rain, precip_snow, snow, qsurf, & |
|---|
| 785 | radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & |
|---|
| 786 | f_qsat_oce, AcoefH, AcoefQ, BcoefH, BcoefQ, & |
|---|
| 787 | tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) |
|---|
| 788 | DO i = 1, knon |
|---|
| 789 | IF (snow(i)<snow_min) tice(knindex(i)) = tsurf_new(i) |
|---|
| 790 | END DO |
|---|
| 791 | |
|---|
| 792 | ! calcul_flux_wind |
|---|
| 793 | CALL calcul_flux_wind(knon, dtime, & |
|---|
| 794 | u0, v0, u1, v1, gustiness, cdragm, & |
|---|
| 795 | AcoefU, AcoefV, BcoefU, BcoefV, & |
|---|
| 796 | p1lay, temp_air, & |
|---|
| 797 | flux_u1, flux_v1) |
|---|
| 798 | |
|---|
| 799 | !**************************************************************************************** |
|---|
| 800 | ! 2) Update snow and ice surface |
|---|
| 801 | !**************************************************************************************** |
|---|
| 802 | ! snow precip |
|---|
| 803 | DO i = 1, knon |
|---|
| 804 | ki = knindex(i) |
|---|
| 805 | IF (precip_snow(i) > 0.) THEN |
|---|
| 806 | snow(i) = snow(i) + precip_snow(i) * dtime * (1. - snow_wfact * (1. - fsic(ki))) |
|---|
| 807 | END IF |
|---|
| 808 | ! snow and ice sublimation |
|---|
| 809 | IF (evap(i) > 0.) THEN |
|---|
| 810 | snow_evap = MIN (snow(i) / dtime, evap(i)) |
|---|
| 811 | snow(i) = snow(i) - snow_evap * dtime |
|---|
| 812 | snow(i) = MAX(0.0, snow(i)) |
|---|
| 813 | seaice(ki) = MAX(0.0, seaice(ki) - (evap(i) - snow_evap) * dtime) |
|---|
| 814 | ENDIF |
|---|
| 815 | ! Melt / Freeze snow from above if Tsurf>0 |
|---|
| 816 | IF (tsurf_new(i)>t_melt) THEN |
|---|
| 817 | ! energy available for melting snow (in kg of melted snow /m2) |
|---|
| 818 | e_melt = MIN(MAX(snow(i) * (tsurf_new(i) - t_melt) * ice_cap / 2. & |
|---|
| 819 | / (ice_lat + ice_cap / 2. * (t_melt - tice(ki))), 0.0), snow(i)) |
|---|
| 820 | ! remove snow |
|---|
| 821 | IF (snow(i)>e_melt) THEN |
|---|
| 822 | snow(i) = snow(i) - e_melt |
|---|
| 823 | tsurf_new(i) = t_melt |
|---|
| 824 | ELSE ! all snow is melted |
|---|
| 825 | ! add remaining heat flux to ice |
|---|
| 826 | e_melt = e_melt - snow(i) |
|---|
| 827 | tice(ki) = tice(ki) + e_melt * ice_lat * 2. / (ice_cap * seaice(ki)) |
|---|
| 828 | tsurf_new(i) = tice(ki) |
|---|
| 829 | END IF |
|---|
| 830 | END IF |
|---|
| 831 | ! melt ice from above if Tice>0 |
|---|
| 832 | IF (tice(ki)>t_melt) THEN |
|---|
| 833 | ! quantity of ice melted (kg/m2) |
|---|
| 834 | e_melt = MAX(seaice(ki) * (tice(ki) - t_melt) * ice_cap / 2. & |
|---|
| 835 | / (ice_lat + ice_cap / 2. * (t_melt - t_freeze)), 0.0) |
|---|
| 836 | ! melt from above, height only |
|---|
| 837 | dhsic = MIN(seaice(ki) - h_ice_min, e_melt) |
|---|
| 838 | e_melt = e_melt - dhsic |
|---|
| 839 | IF (e_melt>0) THEN |
|---|
| 840 | ! lateral melt if ice too thin |
|---|
| 841 | dfsic = MAX(fsic(ki) - ice_frac_min, e_melt / h_ice_min * fsic(ki)) |
|---|
| 842 | ! if all melted add remaining heat to ocean |
|---|
| 843 | e_melt = MAX(0., e_melt * fsic(ki) - dfsic * h_ice_min) |
|---|
| 844 | slab_bilg(ki) = slab_bilg(ki) + e_melt * ice_lat / dtime |
|---|
| 845 | ! update height and fraction |
|---|
| 846 | fsic(ki) = fsic(ki) - dfsic |
|---|
| 847 | END IF |
|---|
| 848 | seaice(ki) = seaice(ki) - dhsic |
|---|
| 849 | ! surface temperature at melting point |
|---|
| 850 | tice(ki) = t_melt |
|---|
| 851 | tsurf_new(i) = t_melt |
|---|
| 852 | END IF |
|---|
| 853 | ! convert snow to ice if below floating line |
|---|
| 854 | h_test = (seaice(ki) + snow(i)) * ice_den - seaice(ki) * sea_den |
|---|
| 855 | IF (h_test>0.) THEN !snow under water |
|---|
| 856 | ! extra snow converted to ice (with added frozen sea water) |
|---|
| 857 | dhsic = h_test / (sea_den - ice_den + sno_den) |
|---|
| 858 | seaice(ki) = seaice(ki) + dhsic |
|---|
| 859 | snow(i) = snow(i) - dhsic * sno_den / ice_den |
|---|
| 860 | ! available energy (freeze sea water + bring to tice) |
|---|
| 861 | e_melt = dhsic * (1. - sno_den / ice_den) * (ice_lat + & |
|---|
| 862 | ice_cap / 2. * (t_freeze - tice(ki))) |
|---|
| 863 | ! update ice temperature |
|---|
| 864 | tice(ki) = tice(ki) + 2. * e_melt / ice_cap / (snow(i) + seaice(ki)) |
|---|
| 865 | END IF |
|---|
| 866 | END DO |
|---|
| 867 | |
|---|
| 868 | ! New albedo |
|---|
| 869 | DO i = 1, knon |
|---|
| 870 | ki = knindex(i) |
|---|
| 871 | ! snow albedo: update snow age |
|---|
| 872 | IF (snow(i)>0.0001) THEN |
|---|
| 873 | agesno(i) = (agesno(i) + (1. - agesno(i) / 50.) * dtime / 86400.)& |
|---|
| 874 | * EXP(-1. * MAX(0.0, precip_snow(i)) * dtime / 5.) |
|---|
| 875 | ELSE |
|---|
| 876 | agesno(i) = 0.0 |
|---|
| 877 | END IF |
|---|
| 878 | ! snow albedo |
|---|
| 879 | alb_snow = alb_sno_min + alb_sno_del * EXP(-agesno(i) / 50.) |
|---|
| 880 | ! ice albedo (varies with ice tkickness and temp) |
|---|
| 881 | alb_ice = MAX(0.0, 0.13 * LOG(100. * seaice(ki) / ice_den) + 0.1) |
|---|
| 882 | IF (tice(ki)>t_freeze - 0.01) THEN |
|---|
| 883 | alb_ice = MIN(alb_ice, alb_ice_wet) |
|---|
| 884 | ELSE |
|---|
| 885 | alb_ice = MIN(alb_ice, alb_ice_dry) |
|---|
| 886 | END IF |
|---|
| 887 | ! pond albedo |
|---|
| 888 | alb_pond = 0.36 - 0.1 * (2.0 + MIN(0.0, MAX(tice(ki) - t_melt, -2.0))) |
|---|
| 889 | ! pond fraction |
|---|
| 890 | frac_pond = 0.2 * (2.0 + MIN(0.0, MAX(tice(ki) - t_melt, -2.0))) |
|---|
| 891 | ! snow fraction |
|---|
| 892 | frac_snow = MAX(0.0, MIN(1.0 - frac_pond, snow(i) / snow_min)) |
|---|
| 893 | ! ice fraction |
|---|
| 894 | frac_ice = MAX(0.0, 1. - frac_pond - frac_snow) |
|---|
| 895 | ! total albedo |
|---|
| 896 | alb1_new(i) = alb_snow * frac_snow + alb_ice * frac_ice + alb_pond * frac_pond |
|---|
| 897 | END DO |
|---|
| 898 | alb2_new(:) = alb1_new(:) |
|---|
| 899 | |
|---|
| 900 | !**************************************************************************************** |
|---|
| 901 | ! 3) Recalculate new ocean temperature (add fluxes below ice) |
|---|
| 902 | ! Melt / freeze from below |
|---|
| 903 | !***********************************************o***************************************** |
|---|
| 904 | !cumul fluxes |
|---|
| 905 | bilg_cum(:) = bilg_cum(:) + slab_bilg(:) |
|---|
| 906 | IF (MOD(itime, cpl_pas)==0) THEN ! time to update tslab & fraction |
|---|
| 907 | ! Add cumulated surface fluxes |
|---|
| 908 | tslab(:, 1) = tslab(:, 1) + bilg_cum(:) * cyang * dtime |
|---|
| 909 | DO i = 1, knon |
|---|
| 910 | ki = knindex(i) |
|---|
| 911 | ! split lateral/top melt-freeze |
|---|
| 912 | frac_mf = MIN(1., MAX(0., (seaice(ki) - h_ice_thin) / (h_ice_thick - h_ice_thin))) |
|---|
| 913 | IF (tslab(ki, 1)<=t_freeze) THEN |
|---|
| 914 | ! ****** Form new ice from below ******* |
|---|
| 915 | ! quantity of new ice |
|---|
| 916 | e_melt = (t_freeze - tslab(ki, 1)) / cyang & |
|---|
| 917 | / (ice_lat + ice_cap / 2. * (t_freeze - tice(ki))) |
|---|
| 918 | ! first increase height to h_thin |
|---|
| 919 | dhsic = MAX(0., MIN(h_ice_thin - seaice(ki), e_melt / fsic(ki))) |
|---|
| 920 | seaice(ki) = dhsic + seaice(ki) |
|---|
| 921 | e_melt = e_melt - fsic(ki) * dhsic |
|---|
| 922 | IF (e_melt>0.) THEN |
|---|
| 923 | ! frac_mf fraction used for lateral increase |
|---|
| 924 | dfsic = MIN(ice_frac_max - fsic(ki), e_melt * frac_mf / seaice(ki)) |
|---|
| 925 | fsic(ki) = fsic(ki) + dfsic |
|---|
| 926 | e_melt = e_melt - dfsic * seaice(ki) |
|---|
| 927 | ! rest used to increase height |
|---|
| 928 | seaice(ki) = MIN(h_ice_max, seaice(ki) + e_melt / fsic(ki)) |
|---|
| 929 | END IF |
|---|
| 930 | tslab(ki, 1) = t_freeze |
|---|
| 931 | ELSE ! slab temperature above freezing |
|---|
| 932 | ! ****** melt ice from below ******* |
|---|
| 933 | ! quantity of melted ice |
|---|
| 934 | e_melt = (tslab(ki, 1) - t_freeze) / cyang & |
|---|
| 935 | / (ice_lat + ice_cap / 2. * (tice(ki) - t_freeze)) |
|---|
| 936 | ! first decrease height to h_thick |
|---|
| 937 | dhsic = MAX(0., MIN(seaice(ki) - h_ice_thick, e_melt / fsic(ki))) |
|---|
| 938 | seaice(ki) = seaice(ki) - dhsic |
|---|
| 939 | e_melt = e_melt - fsic(ki) * dhsic |
|---|
| 940 | IF (e_melt>0) THEN |
|---|
| 941 | ! frac_mf fraction used for height decrease |
|---|
| 942 | dhsic = MAX(0., MIN(seaice(ki) - h_ice_min, e_melt * frac_mf / fsic(ki))) |
|---|
| 943 | seaice(ki) = seaice(ki) - dhsic |
|---|
| 944 | e_melt = e_melt - fsic(ki) * dhsic |
|---|
| 945 | ! rest used to decrease fraction (up to 0!) |
|---|
| 946 | dfsic = MIN(fsic(ki), e_melt / seaice(ki)) |
|---|
| 947 | ! keep remaining in ocean |
|---|
| 948 | e_melt = e_melt - dfsic * seaice(ki) |
|---|
| 949 | END IF |
|---|
| 950 | tslab(ki, 1) = t_freeze + e_melt * ice_lat * cyang |
|---|
| 951 | fsic(ki) = fsic(ki) - dfsic |
|---|
| 952 | END IF |
|---|
| 953 | END DO |
|---|
| 954 | bilg_cum(:) = 0. |
|---|
| 955 | END IF ! coupling time |
|---|
| 956 | |
|---|
| 957 | !tests ice fraction |
|---|
| 958 | WHERE (fsic<ice_frac_min) |
|---|
| 959 | tslab(:, 1) = tslab(:, 1) - fsic * seaice * ice_lat * cyang |
|---|
| 960 | tice = t_melt |
|---|
| 961 | fsic = 0. |
|---|
| 962 | seaice = 0. |
|---|
| 963 | END WHERE |
|---|
| 964 | |
|---|
| 965 | END SUBROUTINE ocean_slab_ice |
|---|
| 966 | |
|---|
| 967 | !**************************************************************************************** |
|---|
| 968 | |
|---|
| 969 | SUBROUTINE ocean_slab_final |
|---|
| 970 | |
|---|
| 971 | !**************************************************************************************** |
|---|
| 972 | ! Deallocate module variables |
|---|
| 973 | !**************************************************************************************** |
|---|
| 974 | IF (ALLOCATED(tslab)) DEALLOCATE(tslab) |
|---|
| 975 | IF (ALLOCATED(fsic)) DEALLOCATE(fsic) |
|---|
| 976 | IF (ALLOCATED(tice)) DEALLOCATE(tice) |
|---|
| 977 | IF (ALLOCATED(seaice)) DEALLOCATE(seaice) |
|---|
| 978 | IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg) |
|---|
| 979 | IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum) |
|---|
| 980 | IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum) |
|---|
| 981 | IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum) |
|---|
| 982 | IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum) |
|---|
| 983 | IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman) |
|---|
| 984 | IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff) |
|---|
| 985 | IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm) |
|---|
| 986 | IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux) |
|---|
| 987 | |
|---|
| 988 | END SUBROUTINE ocean_slab_final |
|---|
| 989 | |
|---|
| 990 | END MODULE ocean_slab_mod |
|---|
