| Line |
Branch |
Exec |
Source |
| 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 |
|
|
|
| 8 |
|
|
USE dimphy |
| 9 |
|
|
USE indice_sol_mod |
| 10 |
|
|
USE surface_data |
| 11 |
|
|
USE mod_grid_phy_lmdz, ONLY: klon_glo |
| 12 |
|
|
USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root |
| 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 ioipsl_getin_p_mod, ONLY : getin_p |
| 132 |
|
|
USE mod_phys_lmdz_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.GT.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.LT.2) THEN |
| 189 |
|
✗ |
slab_gm=.FALSE. |
| 190 |
|
|
ENDIF |
| 191 |
|
|
! Convective adjustment |
| 192 |
|
✗ |
IF (nslay.EQ.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.GT.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.GT.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 mod_phys_lmdz_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.GT.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.GT.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).EQ.0) THEN ! time to update tslab & fraction |
| 490 |
|
|
! *********************************** |
| 491 |
|
|
! Horizontal transport |
| 492 |
|
|
! *********************************** |
| 493 |
|
✗ |
IF (slab_ekman.GT.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.GT.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.GT.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.GT.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.GT.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.GT.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.GT.1).AND.(slab_cadj.GT.0)) THEN |
| 600 |
|
✗ |
DO i=1,klon |
| 601 |
|
✗ |
IF (tslab(i,2).GT.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).LT.t_freeze).OR.(fsic(ki).GT.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).LT.epsfra) THEN ! Free of ice |
| 631 |
|
✗ |
IF (tslab(ki,1).LT.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).GT.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).LT.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 |
|
|
|
| 754 |
|
✗ |
DO i=1,knon |
| 755 |
|
✗ |
ki=knindex(i) |
| 756 |
|
✗ |
IF (snow(i).GT.snow_min) THEN |
| 757 |
|
|
! snow-layer heat capacity |
| 758 |
|
✗ |
cal(i)=2.*RCPD/(snow(i)*ice_cap) |
| 759 |
|
|
! snow conductive flux |
| 760 |
|
✗ |
f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i) |
| 761 |
|
|
! all shortwave flux absorbed |
| 762 |
|
|
f_swpen=0. |
| 763 |
|
|
! bottom flux (ice conduction) |
| 764 |
|
✗ |
slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki) |
| 765 |
|
|
! update ice temperature |
| 766 |
|
|
tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) & |
| 767 |
|
✗ |
*(slab_bilg(ki)+f_cond)*dtime |
| 768 |
|
|
ELSE ! bare ice |
| 769 |
|
|
! ice-layer heat capacity |
| 770 |
|
✗ |
cal(i)=2.*RCPD/(seaice(ki)*ice_cap) |
| 771 |
|
|
! conductive flux |
| 772 |
|
✗ |
f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki) |
| 773 |
|
|
! penetrative shortwave flux... |
| 774 |
|
✗ |
f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den) |
| 775 |
|
✗ |
slab_bilg(ki)=f_swpen-f_cond |
| 776 |
|
|
END IF |
| 777 |
|
✗ |
radsol(i)=radsol(i)+f_cond-f_swpen |
| 778 |
|
|
END DO |
| 779 |
|
|
! weight fluxes to ocean by sea ice fraction |
| 780 |
|
✗ |
slab_bilg(:)=slab_bilg(:)*fsic(:) |
| 781 |
|
|
|
| 782 |
|
|
! calcul_fluxs (sens, lat etc) |
| 783 |
|
|
CALL calcul_fluxs(knon, is_sic, dtime, & |
| 784 |
|
|
tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, & |
| 785 |
|
|
precip_rain, precip_snow, snow, qsurf, & |
| 786 |
|
|
radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & |
| 787 |
|
|
f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, & |
| 788 |
|
✗ |
tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l) |
| 789 |
|
✗ |
DO i=1,knon |
| 790 |
|
✗ |
IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i) |
| 791 |
|
|
END DO |
| 792 |
|
|
|
| 793 |
|
|
! calcul_flux_wind |
| 794 |
|
|
CALL calcul_flux_wind(knon, dtime, & |
| 795 |
|
|
u0, v0, u1, v1, gustiness, cdragm, & |
| 796 |
|
|
AcoefU, AcoefV, BcoefU, BcoefV, & |
| 797 |
|
|
p1lay, temp_air, & |
| 798 |
|
✗ |
flux_u1, flux_v1) |
| 799 |
|
|
|
| 800 |
|
|
!**************************************************************************************** |
| 801 |
|
|
! 2) Update snow and ice surface |
| 802 |
|
|
!**************************************************************************************** |
| 803 |
|
|
! snow precip |
| 804 |
|
✗ |
DO i=1,knon |
| 805 |
|
✗ |
ki=knindex(i) |
| 806 |
|
✗ |
IF (precip_snow(i) > 0.) THEN |
| 807 |
|
✗ |
snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki))) |
| 808 |
|
|
END IF |
| 809 |
|
|
! snow and ice sublimation |
| 810 |
|
✗ |
IF (evap(i) > 0.) THEN |
| 811 |
|
✗ |
snow_evap = MIN (snow(i) / dtime, evap(i)) |
| 812 |
|
✗ |
snow(i) = snow(i) - snow_evap * dtime |
| 813 |
|
✗ |
snow(i) = MAX(0.0, snow(i)) |
| 814 |
|
✗ |
seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime) |
| 815 |
|
|
ENDIF |
| 816 |
|
|
! Melt / Freeze snow from above if Tsurf>0 |
| 817 |
|
✗ |
IF (tsurf_new(i).GT.t_melt) THEN |
| 818 |
|
|
! energy available for melting snow (in kg of melted snow /m2) |
| 819 |
|
|
e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. & |
| 820 |
|
✗ |
/(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i)) |
| 821 |
|
|
! remove snow |
| 822 |
|
✗ |
IF (snow(i).GT.e_melt) THEN |
| 823 |
|
✗ |
snow(i)=snow(i)-e_melt |
| 824 |
|
✗ |
tsurf_new(i)=t_melt |
| 825 |
|
|
ELSE ! all snow is melted |
| 826 |
|
|
! add remaining heat flux to ice |
| 827 |
|
✗ |
e_melt=e_melt-snow(i) |
| 828 |
|
✗ |
tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki)) |
| 829 |
|
✗ |
tsurf_new(i)=tice(ki) |
| 830 |
|
|
END IF |
| 831 |
|
|
END IF |
| 832 |
|
|
! melt ice from above if Tice>0 |
| 833 |
|
✗ |
IF (tice(ki).GT.t_melt) THEN |
| 834 |
|
|
! quantity of ice melted (kg/m2) |
| 835 |
|
|
e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. & |
| 836 |
|
✗ |
/(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0) |
| 837 |
|
|
! melt from above, height only |
| 838 |
|
✗ |
dhsic=MIN(seaice(ki)-h_ice_min,e_melt) |
| 839 |
|
✗ |
e_melt=e_melt-dhsic |
| 840 |
|
✗ |
IF (e_melt.GT.0) THEN |
| 841 |
|
|
! lateral melt if ice too thin |
| 842 |
|
✗ |
dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki)) |
| 843 |
|
|
! if all melted add remaining heat to ocean |
| 844 |
|
✗ |
e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min) |
| 845 |
|
✗ |
slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime |
| 846 |
|
|
! update height and fraction |
| 847 |
|
✗ |
fsic(ki)=fsic(ki)-dfsic |
| 848 |
|
|
END IF |
| 849 |
|
✗ |
seaice(ki)=seaice(ki)-dhsic |
| 850 |
|
|
! surface temperature at melting point |
| 851 |
|
✗ |
tice(ki)=t_melt |
| 852 |
|
✗ |
tsurf_new(i)=t_melt |
| 853 |
|
|
END IF |
| 854 |
|
|
! convert snow to ice if below floating line |
| 855 |
|
✗ |
h_test=(seaice(ki)+snow(i))*ice_den-seaice(ki)*sea_den |
| 856 |
|
✗ |
IF (h_test.GT.0.) THEN !snow under water |
| 857 |
|
|
! extra snow converted to ice (with added frozen sea water) |
| 858 |
|
✗ |
dhsic=h_test/(sea_den-ice_den+sno_den) |
| 859 |
|
✗ |
seaice(ki)=seaice(ki)+dhsic |
| 860 |
|
✗ |
snow(i)=snow(i)-dhsic*sno_den/ice_den |
| 861 |
|
|
! available energy (freeze sea water + bring to tice) |
| 862 |
|
|
e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ & |
| 863 |
|
✗ |
ice_cap/2.*(t_freeze-tice(ki))) |
| 864 |
|
|
! update ice temperature |
| 865 |
|
✗ |
tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki)) |
| 866 |
|
|
END IF |
| 867 |
|
|
END DO |
| 868 |
|
|
|
| 869 |
|
|
! New albedo |
| 870 |
|
✗ |
DO i=1,knon |
| 871 |
|
✗ |
ki=knindex(i) |
| 872 |
|
|
! snow albedo: update snow age |
| 873 |
|
✗ |
IF (snow(i).GT.0.0001) THEN |
| 874 |
|
|
agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)& |
| 875 |
|
✗ |
* EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.) |
| 876 |
|
|
ELSE |
| 877 |
|
✗ |
agesno(i)=0.0 |
| 878 |
|
|
END IF |
| 879 |
|
|
! snow albedo |
| 880 |
|
✗ |
alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.) |
| 881 |
|
|
! ice albedo (varies with ice tkickness and temp) |
| 882 |
|
✗ |
alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1) |
| 883 |
|
✗ |
IF (tice(ki).GT.t_freeze-0.01) THEN |
| 884 |
|
✗ |
alb_ice=MIN(alb_ice,alb_ice_wet) |
| 885 |
|
|
ELSE |
| 886 |
|
✗ |
alb_ice=MIN(alb_ice,alb_ice_dry) |
| 887 |
|
|
END IF |
| 888 |
|
|
! pond albedo |
| 889 |
|
✗ |
alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) |
| 890 |
|
|
! pond fraction |
| 891 |
|
✗ |
frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0))) |
| 892 |
|
|
! snow fraction |
| 893 |
|
✗ |
frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min)) |
| 894 |
|
|
! ice fraction |
| 895 |
|
✗ |
frac_ice=MAX(0.0,1.-frac_pond-frac_snow) |
| 896 |
|
|
! total albedo |
| 897 |
|
✗ |
alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond |
| 898 |
|
|
END DO |
| 899 |
|
✗ |
alb2_new(:) = alb1_new(:) |
| 900 |
|
|
|
| 901 |
|
|
!**************************************************************************************** |
| 902 |
|
|
! 3) Recalculate new ocean temperature (add fluxes below ice) |
| 903 |
|
|
! Melt / freeze from below |
| 904 |
|
|
!***********************************************o***************************************** |
| 905 |
|
|
!cumul fluxes |
| 906 |
|
✗ |
bilg_cum(:)=bilg_cum(:)+slab_bilg(:) |
| 907 |
|
✗ |
IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction |
| 908 |
|
|
! Add cumulated surface fluxes |
| 909 |
|
✗ |
tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime |
| 910 |
|
✗ |
DO i=1,knon |
| 911 |
|
✗ |
ki=knindex(i) |
| 912 |
|
|
! split lateral/top melt-freeze |
| 913 |
|
✗ |
frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin))) |
| 914 |
|
✗ |
IF (tslab(ki,1).LE.t_freeze) THEN |
| 915 |
|
|
! ****** Form new ice from below ******* |
| 916 |
|
|
! quantity of new ice |
| 917 |
|
|
e_melt=(t_freeze-tslab(ki,1))/cyang & |
| 918 |
|
✗ |
/(ice_lat+ice_cap/2.*(t_freeze-tice(ki))) |
| 919 |
|
|
! first increase height to h_thin |
| 920 |
|
✗ |
dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki))) |
| 921 |
|
✗ |
seaice(ki)=dhsic+seaice(ki) |
| 922 |
|
✗ |
e_melt=e_melt-fsic(ki)*dhsic |
| 923 |
|
✗ |
IF (e_melt.GT.0.) THEN |
| 924 |
|
|
! frac_mf fraction used for lateral increase |
| 925 |
|
✗ |
dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki)) |
| 926 |
|
✗ |
fsic(ki)=fsic(ki)+dfsic |
| 927 |
|
✗ |
e_melt=e_melt-dfsic*seaice(ki) |
| 928 |
|
|
! rest used to increase height |
| 929 |
|
✗ |
seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki)) |
| 930 |
|
|
END IF |
| 931 |
|
✗ |
tslab(ki,1)=t_freeze |
| 932 |
|
|
ELSE ! slab temperature above freezing |
| 933 |
|
|
! ****** melt ice from below ******* |
| 934 |
|
|
! quantity of melted ice |
| 935 |
|
|
e_melt=(tslab(ki,1)-t_freeze)/cyang & |
| 936 |
|
✗ |
/(ice_lat+ice_cap/2.*(tice(ki)-t_freeze)) |
| 937 |
|
|
! first decrease height to h_thick |
| 938 |
|
✗ |
dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki))) |
| 939 |
|
✗ |
seaice(ki)=seaice(ki)-dhsic |
| 940 |
|
✗ |
e_melt=e_melt-fsic(ki)*dhsic |
| 941 |
|
✗ |
IF (e_melt.GT.0) THEN |
| 942 |
|
|
! frac_mf fraction used for height decrease |
| 943 |
|
✗ |
dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki))) |
| 944 |
|
✗ |
seaice(ki)=seaice(ki)-dhsic |
| 945 |
|
✗ |
e_melt=e_melt-fsic(ki)*dhsic |
| 946 |
|
|
! rest used to decrease fraction (up to 0!) |
| 947 |
|
✗ |
dfsic=MIN(fsic(ki),e_melt/seaice(ki)) |
| 948 |
|
|
! keep remaining in ocean |
| 949 |
|
✗ |
e_melt=e_melt-dfsic*seaice(ki) |
| 950 |
|
|
END IF |
| 951 |
|
✗ |
tslab(ki,1)=t_freeze+e_melt*ice_lat*cyang |
| 952 |
|
✗ |
fsic(ki)=fsic(ki)-dfsic |
| 953 |
|
|
END IF |
| 954 |
|
|
END DO |
| 955 |
|
✗ |
bilg_cum(:)=0. |
| 956 |
|
|
END IF ! coupling time |
| 957 |
|
|
|
| 958 |
|
|
!tests ice fraction |
| 959 |
|
✗ |
WHERE (fsic.LT.ice_frac_min) |
| 960 |
|
|
tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang |
| 961 |
|
|
tice=t_melt |
| 962 |
|
|
fsic=0. |
| 963 |
|
|
seaice=0. |
| 964 |
|
|
END WHERE |
| 965 |
|
|
|
| 966 |
|
✗ |
END SUBROUTINE ocean_slab_ice |
| 967 |
|
|
! |
| 968 |
|
|
!**************************************************************************************** |
| 969 |
|
|
! |
| 970 |
|
|
SUBROUTINE ocean_slab_final |
| 971 |
|
|
|
| 972 |
|
|
!**************************************************************************************** |
| 973 |
|
|
! Deallocate module variables |
| 974 |
|
|
!**************************************************************************************** |
| 975 |
|
|
IF (ALLOCATED(tslab)) DEALLOCATE(tslab) |
| 976 |
|
|
IF (ALLOCATED(fsic)) DEALLOCATE(fsic) |
| 977 |
|
|
IF (ALLOCATED(tice)) DEALLOCATE(tice) |
| 978 |
|
|
IF (ALLOCATED(seaice)) DEALLOCATE(seaice) |
| 979 |
|
|
IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg) |
| 980 |
|
|
IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum) |
| 981 |
|
|
IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum) |
| 982 |
|
|
IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum) |
| 983 |
|
|
IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum) |
| 984 |
|
|
IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman) |
| 985 |
|
|
IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff) |
| 986 |
|
|
IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm) |
| 987 |
|
|
IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux) |
| 988 |
|
|
|
| 989 |
|
✗ |
END SUBROUTINE ocean_slab_final |
| 990 |
|
|
|
| 991 |
|
|
END MODULE ocean_slab_mod |
| 992 |
|
|
|