LMDZ
ocean_slab_mod.F90
Go to the documentation of this file.
1 !
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 
12  IMPLICIT NONE
13  PRIVATE
14  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
15 
16 !****************************************************************************************
17 ! Global saved variables
18 !****************************************************************************************
19 
20  INTEGER, PRIVATE, SAVE :: cpl_pas
21  !$OMP THREADPRIVATE(cpl_pas)
22  REAL, PRIVATE, SAVE :: cyang
23  !$OMP THREADPRIVATE(cyang)
24  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: slabh
25  !$OMP THREADPRIVATE(slabh)
26  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE :: tslab
27  !$OMP THREADPRIVATE(tslab)
28  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: fsic
29  !$OMP THREADPRIVATE(fsic)
30  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: tice
31  !$OMP THREADPRIVATE(tice)
32  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: seaice
33  !$OMP THREADPRIVATE(seaice)
34  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bils
35  !$OMP THREADPRIVATE(slab_bils)
36  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bils_cum
37  !$OMP THREADPRIVATE(bils_cum)
38  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE :: slab_bilg
39  !$OMP THREADPRIVATE(slab_bilg)
40  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: bilg_cum
41  !$OMP THREADPRIVATE(bilg_cum)
42 
43 !****************************************************************************************
44 ! Parameters (could be read in def file: move to slab_init)
45 !****************************************************************************************
46 ! snow and ice physical characteristics:
47  REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
48  REAL, PARAMETER :: t_melt=273.15 ! melting ice temp
49  REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
50  REAL, PARAMETER :: ice_den=917.
51  REAL, PARAMETER :: sea_den=1025.
52  REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
53  REAL, PARAMETER :: sno_cond=0.31*sno_den
54  REAL, PARAMETER :: ice_cap=2067. ! specific heat capacity, snow and ice
55  REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
56 
57 ! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
58  REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm
59  REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean
60  REAL, PARAMETER :: ice_frac_min=0.001
61  REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction
62  REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness
63  REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness
64  ! below ice_thin, priority is melt lateral / grow height
65  ! ice_thin is also height of new ice
66  REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness
67  ! above ice_thick, priority is melt height / grow lateral
68  REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice
69  REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height
70 
71 ! albedo and radiation parameters
72  REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo
73  REAL, PARAMETER :: alb_sno_del=0.3 !max snow albedo = min + del
74  REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
75  REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
76  REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
77  REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
78 
79 !****************************************************************************************
80 
81 CONTAINS
82 !
83 !****************************************************************************************
84 !
85  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
86  !, seaice_rst etc
87 
88  use ioipsl
89 
90  ! For ok_xxx vars (Ekman...)
91  include "clesphys.h"
92 
93  ! Input variables
94 !****************************************************************************************
95  REAL, INTENT(IN) :: dtime
96 ! Variables read from restart file
97  REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
98 
99 ! Local variables
100 !****************************************************************************************
101  INTEGER :: error
102  CHARACTER (len = 80) :: abort_message
103  CHARACTER (len = 20) :: modname = 'ocean_slab_intit'
104 
105 !****************************************************************************************
106 ! Allocate surface fraction read from restart file
107 !****************************************************************************************
108  ALLOCATE(fsic(klon), stat = error)
109  IF (error /= 0) THEN
110  abort_message='Pb allocation tmp_pctsrf_slab'
111  CALL abort_physic(modname,abort_message,1)
112  ENDIF
113  fsic(:)=0.
114  WHERE (1.-zmasq(:)>epsfra)
115  fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
116  END WHERE
117 
118 !****************************************************************************************
119 ! Allocate saved variables
120 !****************************************************************************************
121  ALLOCATE(tslab(klon,nslay), stat=error)
122  IF (error /= 0) CALL abort_physic &
123  (modname,'pb allocation tslab', 1)
124 
125  ALLOCATE(slab_bils(klon), stat = error)
126  IF (error /= 0) THEN
127  abort_message='Pb allocation slab_bils'
128  CALL abort_physic(modname,abort_message,1)
129  ENDIF
130  slab_bils(:) = 0.0
131  ALLOCATE(bils_cum(klon), stat = error)
132  IF (error /= 0) THEN
133  abort_message='Pb allocation slab_bils_cum'
134  CALL abort_physic(modname,abort_message,1)
135  ENDIF
136  bils_cum(:) = 0.0
137 
138  IF (version_ocean=='sicINT') THEN
139  ALLOCATE(slab_bilg(klon), stat = error)
140  IF (error /= 0) THEN
141  abort_message='Pb allocation slab_bilg'
142  CALL abort_physic(modname,abort_message,1)
143  ENDIF
144  slab_bilg(:) = 0.0
145  ALLOCATE(bilg_cum(klon), stat = error)
146  IF (error /= 0) THEN
147  abort_message='Pb allocation slab_bilg_cum'
148  CALL abort_physic(modname,abort_message,1)
149  ENDIF
150  bilg_cum(:) = 0.0
151  ALLOCATE(tice(klon), stat = error)
152  IF (error /= 0) THEN
153  abort_message='Pb allocation slab_tice'
154  CALL abort_physic(modname,abort_message,1)
155  ENDIF
156  ALLOCATE(seaice(klon), stat = error)
157  IF (error /= 0) THEN
158  abort_message='Pb allocation slab_seaice'
159  CALL abort_physic(modname,abort_message,1)
160  ENDIF
161  END IF
162 
163 !****************************************************************************************
164 ! Define some parameters
165 !****************************************************************************************
166 ! Layer thickness
167  ALLOCATE(slabh(nslay), stat = error)
168  IF (error /= 0) THEN
169  abort_message='Pb allocation slabh'
170  CALL abort_physic(modname,abort_message,1)
171  ENDIF
172  slabh(1)=50.
173 ! cyang = 1/heat capacity of top layer (rho.c.H)
174  cyang=1/(slabh(1)*4.228e+06)
175 
176 ! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
177 ! pour un calcul à chaque pas de temps, cpl_pas=1
178  cpl_pas = nint(86400./dtime * 1.0) ! une fois par jour
179  CALL getin('cpl_pas',cpl_pas)
180  print *,'cpl_pas',cpl_pas
181 
182  END SUBROUTINE ocean_slab_init
183 !
184 !****************************************************************************************
185 !
186  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
187 
188  USE limit_read_mod
189 
190 ! INCLUDE "clesphys.h"
191 
192 ! Arguments
193 !****************************************************************************************
194  INTEGER, INTENT(IN) :: itime ! numero du pas de temps courant
195  INTEGER, INTENT(IN) :: jour ! jour a lire dans l'annee
196  REAL , INTENT(IN) :: dtime ! pas de temps de la physique (en s)
197  REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg ! sub-surface fraction
198  LOGICAL, INTENT(OUT) :: is_modified ! true if pctsrf is modified at this time step
199 
200 ! Local variables
201 !****************************************************************************************
202 
203  IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN
204  CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
205  ELSE
206  pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
207  pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
208  is_modified=.true.
209  END IF
210 
211  END SUBROUTINE ocean_slab_frac
212 !
213 !****************************************************************************************
214 !
215  SUBROUTINE ocean_slab_noice( &
216  itime, dtime, jour, knon, knindex, &
217  p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, &
218  acoefh, acoefq, bcoefh, bcoefq, &
219  acoefu, acoefv, bcoefu, bcoefv, &
220  ps, u1, v1, gustiness, tsurf_in, &
221  radsol, snow, &
222  qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
223  tsurf_new, dflux_s, dflux_l, qflux)
224 
225  USE calcul_fluxs_mod
226 
227  include "clesphys.h"
228 
229 ! Input arguments
230 !****************************************************************************************
231  INTEGER, INTENT(IN) :: itime
232  INTEGER, INTENT(IN) :: jour
233  INTEGER, INTENT(IN) :: knon
234  INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
235  REAL, INTENT(IN) :: dtime
236  REAL, DIMENSION(klon), INTENT(IN) :: p1lay
237  REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragq, cdragm
238  REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
239  REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
240  REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
241  REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
242  REAL, DIMENSION(klon), INTENT(IN) :: ps
243  REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
244  REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
245  REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
246 
247 ! In/Output arguments
248 !****************************************************************************************
249  REAL, DIMENSION(klon), INTENT(INOUT) :: snow
250 
251 ! Output arguments
252 !****************************************************************************************
253  REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
254  REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
255  REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
256  REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
257  REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
258  REAL, DIMENSION(klon), INTENT(OUT) :: qflux
259 
260 ! Local variables
261 !****************************************************************************************
262  INTEGER :: i,ki
263  ! surface fluxes
264  REAL, DIMENSION(klon) :: cal, beta, dif_grnd
265  ! flux correction
266  REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
267  ! surface wind stress
268  REAL, DIMENSION(klon) :: u0, v0
269  REAL, DIMENSION(klon) :: u1_lay, v1_lay
270  ! ice creation
271  REAL :: e_freeze, h_new, dfsic
272 
273 !****************************************************************************************
274 ! 1) Flux calculation
275 !
276 !****************************************************************************************
277  cal(:) = 0.
278  beta(:) = 1.
279  dif_grnd(:) = 0.
280 
281 ! Suppose zero surface speed
282  u0(:)=0.0
283  v0(:)=0.0
284  u1_lay(:) = u1(:) - u0(:)
285  v1_lay(:) = v1(:) - v0(:)
286 
287  CALL calcul_fluxs(knon, is_oce, dtime, &
288  tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
289  precip_rain, precip_snow, snow, qsurf, &
290  radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
291  f_qsat_oce,acoefh, acoefq, bcoefh, bcoefq, &
292  tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
293 
294 ! - Flux calculation at first modele level for U and V
295  CALL calcul_flux_wind(knon, dtime, &
296  u0, v0, u1, v1, gustiness, cdragm, &
297  acoefu, acoefv, bcoefu, bcoefv, &
298  p1lay, temp_air, &
299  flux_u1, flux_v1)
300 
301 ! Accumulate total fluxes locally
302  slab_bils(:)=0.
303  DO i=1,knon
304  ki=knindex(i)
305  slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
306  -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
307  bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
308 ! Also taux, tauy, saved vars...
309  END DO
310 
311 !****************************************************************************************
312 ! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
313 !
314 !****************************************************************************************
315  lmt_bils(:)=0.
316  CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
317  ! lmt_bils and diff_sst,siv saved by limit_slab
318  qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
319  ! qflux = total QFlux correction (in W/m2)
320 
321 !****************************************************************************************
322 ! 3) Recalculate new temperature
323 !
324 !***********************************************o*****************************************
325  tsurf_new=tsurf_in
326  IF (mod(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
327  ! Compute transport
328  ! Add read QFlux and SST tendency
329  tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
330  ! Add cumulated surface fluxes
331  tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
332  ! Update surface temperature
333  SELECT CASE(version_ocean)
334  CASE('sicNO')
335  DO i=1,knon
336  ki=knindex(i)
337  tsurf_new(i)=tslab(ki,1)
338  END DO
339  CASE('sicOBS') ! check for sea ice or tslab below freezing
340  DO i=1,knon
341  ki=knindex(i)
342  IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
343  tslab(ki,1)=t_freeze
344  END IF
345  tsurf_new(i)=tslab(ki,1)
346  END DO
347  CASE('sicINT')
348  DO i=1,knon
349  ki=knindex(i)
350  IF (fsic(ki).LT.epsfra) THEN ! Free of ice
351  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
352  ! quantity of new ice formed
353  e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
354  ! new ice
355  tice(ki)=t_freeze
356  fsic(ki)=min(ice_frac_max,e_freeze/h_ice_thin)
357  IF (fsic(ki).GT.ice_frac_min) THEN
358  seaice(ki)=min(e_freeze/fsic(ki),h_ice_max)
359  tslab(ki,1)=t_freeze
360  ELSE
361  fsic(ki)=0.
362  END IF
363  tsurf_new(i)=t_freeze
364  ELSE
365  tsurf_new(i)=tslab(ki,1)
366  END IF
367  ELSE ! ice present
368  tsurf_new(i)=t_freeze
369  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
370  ! quantity of new ice formed over open ocean
371  e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
372  /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
373  ! new ice height and fraction
374  h_new=min(h_ice_new,seaice(ki)) ! max new height ice_new
375  dfsic=min(ice_frac_max-fsic(ki),e_freeze/h_new)
376  h_new=min(e_freeze/dfsic,h_ice_max)
377  ! update tslab to freezing over open ocean only
378  tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
379  ! update sea ice
380  seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
381  /(dfsic+fsic(ki))
382  fsic(ki)=fsic(ki)+dfsic
383  ! update snow?
384  END IF !freezing
385  END IF ! sea ice present
386  END DO
387  END SELECT
388  bils_cum(:)=0.0! clear cumulated fluxes
389  END IF ! coupling time
390  END SUBROUTINE ocean_slab_noice
391 !
392 !****************************************************************************************
393 
394  SUBROUTINE ocean_slab_ice( &
395  itime, dtime, jour, knon, knindex, &
396  tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
397  acoefh, acoefq, bcoefh, bcoefq, &
398  acoefu, acoefv, bcoefu, bcoefv, &
399  ps, u1, v1, gustiness, &
400  radsol, snow, qsurf, qsol, agesno, &
401  alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
402  tsurf_new, dflux_s, dflux_l, swnet)
403 
404  USE calcul_fluxs_mod
405 
406  include "YOMCST.h"
407  include "clesphys.h"
408 
409 ! Input arguments
410 !****************************************************************************************
411  INTEGER, INTENT(IN) :: itime, jour, knon
412  INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
413  REAL, INTENT(IN) :: dtime
414  REAL, DIMENSION(klon), INTENT(IN) :: tsurf_in
415  REAL, DIMENSION(klon), INTENT(IN) :: p1lay
416  REAL, DIMENSION(klon), INTENT(IN) :: cdragh, cdragm
417  REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
418  REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum
419  REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ
420  REAL, DIMENSION(klon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV
421  REAL, DIMENSION(klon), INTENT(IN) :: ps
422  REAL, DIMENSION(klon), INTENT(IN) :: u1, v1, gustiness
423  REAL, DIMENSION(klon), INTENT(IN) :: swnet
424 
425 ! In/Output arguments
426 !****************************************************************************************
427  REAL, DIMENSION(klon), INTENT(INOUT) :: snow, qsol
428  REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
429  REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
430 
431 ! Output arguments
432 !****************************************************************************************
433  REAL, DIMENSION(klon), INTENT(OUT) :: qsurf
434  REAL, DIMENSION(klon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval
435  REAL, DIMENSION(klon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval
436  REAL, DIMENSION(klon), INTENT(OUT) :: evap, fluxsens, fluxlat
437  REAL, DIMENSION(klon), INTENT(OUT) :: flux_u1, flux_v1
438  REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new
439  REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s, dflux_l
440 
441 ! Local variables
442 !****************************************************************************************
443  INTEGER :: i,ki
444  REAL, DIMENSION(klon) :: cal, beta, dif_grnd
445  REAL, DIMENSION(klon) :: u0, v0
446  REAL, DIMENSION(klon) :: u1_lay, v1_lay
447  ! intermediate heat fluxes:
448  REAL :: f_cond, f_swpen
449  ! for snow/ice albedo:
450  REAL :: alb_snow, alb_ice, alb_pond
451  REAL :: frac_snow, frac_ice, frac_pond
452  ! for ice melt / freeze
453  REAL :: e_melt, snow_evap, h_test
454  ! dhsic, dfsic change in ice mass, fraction.
455  REAL :: dhsic, dfsic, frac_mf
456 
457 !****************************************************************************************
458 ! 1) Flux calculation
459 !****************************************************************************************
460 ! Suppose zero surface speed
461  u0(:)=0.0
462  v0(:)=0.0
463  u1_lay(:) = u1(:) - u0(:)
464  v1_lay(:) = v1(:) - v0(:)
465 
466 ! set beta, cal, compute conduction fluxes inside ice/snow
467  slab_bilg(:)=0.
468  dif_grnd(:)=0.
469  beta(:) = 1.
470  DO i=1,knon
471  ki=knindex(i)
472  IF (snow(i).GT.snow_min) THEN
473  ! snow-layer heat capacity
474  cal(i)=2.*rcpd/(snow(i)*ice_cap)
475  ! snow conductive flux
476  f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
477  ! all shortwave flux absorbed
478  f_swpen=0.
479  ! bottom flux (ice conduction)
480  slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
481  ! update ice temperature
482  tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
483  *(slab_bilg(ki)+f_cond)*dtime
484  ELSE ! bare ice
485  ! ice-layer heat capacity
486  cal(i)=2.*rcpd/(seaice(ki)*ice_cap)
487  ! conductive flux
488  f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
489  ! penetrative shortwave flux...
490  f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den)
491  slab_bilg(ki)=f_swpen-f_cond
492  END IF
493  radsol(i)=radsol(i)+f_cond-f_swpen
494  END DO
495  ! weight fluxes to ocean by sea ice fraction
496  slab_bilg(:)=slab_bilg(:)*fsic(:)
497 
498 ! calcul_fluxs (sens, lat etc)
499  CALL calcul_fluxs(knon, is_sic, dtime, &
500  tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
501  precip_rain, precip_snow, snow, qsurf, &
502  radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
503  f_qsat_oce,acoefh, acoefq, bcoefh, bcoefq, &
504  tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
505  DO i=1,knon
506  IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
507  END DO
508 
509 ! calcul_flux_wind
510  CALL calcul_flux_wind(knon, dtime, &
511  u0, v0, u1, v1, gustiness, cdragm, &
512  acoefu, acoefv, bcoefu, bcoefv, &
513  p1lay, temp_air, &
514  flux_u1, flux_v1)
515 
516 !****************************************************************************************
517 ! 2) Update snow and ice surface
518 !****************************************************************************************
519 ! snow precip
520  DO i=1,knon
521  ki=knindex(i)
522  IF (precip_snow(i) > 0.) THEN
523  snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
524  END IF
525 ! snow and ice sublimation
526  IF (evap(i) > 0.) THEN
527  snow_evap = min(snow(i) / dtime, evap(i))
528  snow(i) = snow(i) - snow_evap * dtime
529  snow(i) = max(0.0, snow(i))
530  seaice(ki) = max(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
531  ENDIF
532 ! Melt / Freeze from above if Tsurf>0
533  IF (tsurf_new(i).GT.t_melt) THEN
534  ! energy available for melting snow (in kg/m2 of snow)
535  e_melt = min(max(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
536  /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
537  ! remove snow
538  IF (snow(i).GT.e_melt) THEN
539  snow(i)=snow(i)-e_melt
540  tsurf_new(i)=t_melt
541  ELSE ! all snow is melted
542  ! add remaining heat flux to ice
543  e_melt=e_melt-snow(i)
544  tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
545  tsurf_new(i)=tice(ki)
546  END IF
547  END IF
548 ! melt ice from above if Tice>0
549  IF (tice(ki).GT.t_melt) THEN
550  ! quantity of ice melted (kg/m2)
551  e_melt=max(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. &
552  /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
553  ! melt from above, height only
554  dhsic=min(seaice(ki)-h_ice_min,e_melt)
555  e_melt=e_melt-dhsic
556  IF (e_melt.GT.0) THEN
557  ! lateral melt if ice too thin
558  dfsic=max(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
559  ! if all melted add remaining heat to ocean
560  e_melt=max(0.,e_melt*fsic(ki)-dfsic*h_ice_min)
561  slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
562  ! update height and fraction
563  fsic(ki)=fsic(ki)-dfsic
564  END IF
565  seaice(ki)=seaice(ki)-dhsic
566  ! surface temperature at melting point
567  tice(ki)=t_melt
568  tsurf_new(i)=t_melt
569  END IF
570  ! convert snow to ice if below floating line
571  h_test=(seaice(ki)+snow(i))*ice_den-seaice(ki)*sea_den
572  IF (h_test.GT.0.) THEN !snow under water
573  ! extra snow converted to ice (with added frozen sea water)
574  dhsic=h_test/(sea_den-ice_den+sno_den)
575  seaice(ki)=seaice(ki)+dhsic
576  snow(i)=snow(i)-dhsic*sno_den/ice_den
577  ! available energy (freeze sea water + bring to tice)
578  e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
579  ice_cap/2.*(t_freeze-tice(ki)))
580  ! update ice temperature
581  tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
582  END IF
583  END DO
584 
585 ! New albedo
586  DO i=1,knon
587  ki=knindex(i)
588  ! snow albedo: update snow age
589  IF (snow(i).GT.0.0001) THEN
590  agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
591  * exp(-1.*max(0.0,precip_snow(i))*dtime/5.)
592  ELSE
593  agesno(i)=0.0
594  END IF
595  ! snow albedo
596  alb_snow=alb_sno_min+alb_sno_del*exp(-agesno(i)/50.)
597  ! ice albedo (varies with ice tkickness and temp)
598  alb_ice=max(0.0,0.13*log(100.*seaice(ki)/ice_den)+0.1)
599  IF (tice(ki).GT.t_freeze-0.01) THEN
600  alb_ice=min(alb_ice,alb_ice_wet)
601  ELSE
602  alb_ice=min(alb_ice,alb_ice_dry)
603  END IF
604  ! pond albedo
605  alb_pond=0.36-0.1*(2.0+min(0.0,max(tice(ki)-t_melt,-2.0)))
606  ! pond fraction
607  frac_pond=0.2*(2.0+min(0.0,max(tice(ki)-t_melt,-2.0)))
608  ! snow fraction
609  frac_snow=max(0.0,min(1.0-frac_pond,snow(i)/snow_min))
610  ! ice fraction
611  frac_ice=max(0.0,1.-frac_pond-frac_snow)
612  ! total albedo
613  alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond
614  END DO
615  alb2_new(:) = alb1_new(:)
616 
617 !****************************************************************************************
618 ! 3) Recalculate new ocean temperature (add fluxes below ice)
619 ! Melt / freeze from below
620 !***********************************************o*****************************************
621  !cumul fluxes
622  bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
623  IF (mod(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
624  ! Add cumulated surface fluxes
625  tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime
626  DO i=1,knon
627  ki=knindex(i)
628  ! split lateral/top melt-freeze
629  frac_mf=min(1.,max(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
630  IF (tslab(ki,1).LE.t_freeze) THEN
631  ! ****** Form new ice from below *******
632  ! quantity of new ice
633  e_melt=(t_freeze-tslab(ki,1))/cyang &
634  /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
635  ! first increase height to h_thin
636  dhsic=max(0.,min(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
637  seaice(ki)=dhsic+seaice(ki)
638  e_melt=e_melt-fsic(ki)*dhsic
639  IF (e_melt.GT.0.) THEN
640  ! frac_mf fraction used for lateral increase
641  dfsic=min(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
642  fsic(ki)=fsic(ki)+dfsic
643  e_melt=e_melt-dfsic*seaice(ki)
644  ! rest used to increase height
645  seaice(ki)=min(h_ice_max,seaice(ki)+e_melt/fsic(ki))
646  END IF
647  tslab(ki,1)=t_freeze
648  ELSE ! slab temperature above freezing
649  ! ****** melt ice from below *******
650  ! quantity of melted ice
651  e_melt=(tslab(ki,1)-t_freeze)/cyang &
652  /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
653  ! first decrease height to h_thick
654  dhsic=max(0.,min(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
655  seaice(ki)=seaice(ki)-dhsic
656  e_melt=e_melt-fsic(ki)*dhsic
657  IF (e_melt.GT.0) THEN
658  ! frac_mf fraction used for height decrease
659  dhsic=max(0.,min(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
660  seaice(ki)=seaice(ki)-dhsic
661  e_melt=e_melt-fsic(ki)*dhsic
662  ! rest used to decrease fraction (up to 0!)
663  dfsic=min(fsic(ki),e_melt/seaice(ki))
664  ! keep remaining in ocean
665  e_melt=e_melt-dfsic*seaice(ki)
666  END IF
667  tslab(ki,1)=t_freeze+e_melt*ice_lat*cyang
668  fsic(ki)=fsic(ki)-dfsic
669  END IF
670  END DO
671  bilg_cum(:)=0.
672  END IF ! coupling time
673 
674  !tests fraction glace
675  WHERE (fsic.LT.ice_frac_min)
676  tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
677  tice=t_melt
678  fsic=0.
679  seaice=0.
680  END WHERE
681 
682  END SUBROUTINE ocean_slab_ice
683 !
684 !****************************************************************************************
685 !
686  SUBROUTINE ocean_slab_final
687 
688 !****************************************************************************************
689 ! Deallocate module variables
690 !****************************************************************************************
691  IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
692  IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
693  IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
694  IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
695  IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
696  IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
697  IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
698 
699  END SUBROUTINE ocean_slab_final
700 
701 END MODULE ocean_slab_mod
ocean_slab_mod::tice
real, dimension(:), allocatable, save, public tice
Definition: ocean_slab_mod.F90:30
ocean_slab_mod::h_ice_thick
real, parameter h_ice_thick
Definition: ocean_slab_mod.F90:66
limit_read_mod::limit_read_frac
subroutine limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified)
Definition: limit_read_mod.F90:34
ocean_slab_mod::slab_bilg
real, dimension(:), allocatable, save, public slab_bilg
Definition: ocean_slab_mod.F90:38
ocean_slab_mod::ocean_slab_final
subroutine ocean_slab_final
Definition: ocean_slab_mod.F90:687
ocean_slab_mod::t_freeze
real, parameter t_freeze
Definition: ocean_slab_mod.F90:47
ocean_slab_mod::cpl_pas
integer, save, private cpl_pas
Definition: ocean_slab_mod.F90:20
limit_read_mod
Definition: limit_read_mod.F90:4
ocean_slab_mod::ice_cap
real, parameter ice_cap
Definition: ocean_slab_mod.F90:54
ocean_slab_mod::alb_sno_min
real, parameter alb_sno_min
Definition: ocean_slab_mod.F90:72
ocean_slab_mod::h_ice_new
real, parameter h_ice_new
Definition: ocean_slab_mod.F90:68
calcul_fluxs_mod::calcul_flux_wind
subroutine calcul_flux_wind(knon, dtime, u0, v0, u1, v1, gustiness, cdrag_m, AcoefU, AcoefV, BcoefU, BcoefV, p1lay, t1lay, flux_u1, flux_v1)
Definition: calcul_fluxs_mod.F90:258
ocean_slab_mod::t_melt
real, parameter t_melt
Definition: ocean_slab_mod.F90:48
dtime
!$Header!integer nvarmx dtime
Definition: gradsdef.h:20
surface_data::version_ocean
character(len=6), save version_ocean
Definition: surface_data.F90:23
dimphy::klon
integer, save klon
Definition: dimphy.F90:3
ocean_slab_mod::ocean_slab_frac
subroutine, public ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
Definition: ocean_slab_mod.F90:187
ocean_slab_mod::bils_cum
real, dimension(:), allocatable, save, private bils_cum
Definition: ocean_slab_mod.F90:36
ocean_slab_mod::fsic
real, dimension(:), allocatable, save, public fsic
Definition: ocean_slab_mod.F90:28
ocean_slab_mod
Definition: ocean_slab_mod.F90:2
ocean_slab_mod::slabh
real, dimension(:), allocatable, save, private slabh
Definition: ocean_slab_mod.F90:24
f_qsat_oce
!$Id ok_orolf LOGICAL ok_limitvrai LOGICAL ok_all_xml INTEGER iflag_ener_conserv REAL solaire RCFC12 RCFC12_act CFC12_ppt!IM ajout CFMIP2 CMIP5 LOGICAL ok_4xCO2atm RCFC12_per CFC12_ppt_per!OM correction du bilan d eau global!OM Correction sur precip KE REAL cvl_corr!OM Fonte calotte dans bilan eau LOGICAL ok_lic_melt!IM simulateur ISCCP INTEGER overlap!IM seuils cdrh REAL cdhmax!IM param stabilite s terres et en dehors REAL f_ri_cd_min!IM MAFo pmagic evap0!Frottement au f_cdrag_oce REAL f_qsat_oce
Definition: clesphys.h:46
ocean_slab_mod::ice_frac_max
real, parameter ice_frac_max
Definition: ocean_slab_mod.F90:61
limit_slab
subroutine limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
Definition: limit_slab.F90:4
ocean_slab_mod::sno_cond
real, parameter sno_cond
Definition: ocean_slab_mod.F90:53
ocean_slab_mod::pen_frac
real, parameter pen_frac
Definition: ocean_slab_mod.F90:76
ocean_slab_mod::ice_frac_min
real, parameter ice_frac_min
Definition: ocean_slab_mod.F90:60
ocean_slab_mod::alb_sno_del
real, parameter alb_sno_del
Definition: ocean_slab_mod.F90:73
ocean_slab_mod::ice_lat
real, parameter ice_lat
Definition: ocean_slab_mod.F90:55
ocean_slab_mod::ocean_slab_init
subroutine, public ocean_slab_init(dtime, pctsrf_rst)
Definition: ocean_slab_mod.F90:86
dimphy::zmasq
real, dimension(:), allocatable, save zmasq
Definition: dimphy.F90:14
ocean_slab_mod::seaice
real, dimension(:), allocatable, save, public seaice
Definition: ocean_slab_mod.F90:32
ocean_slab_mod::bilg_cum
real, dimension(:), allocatable, save, private bilg_cum
Definition: ocean_slab_mod.F90:40
ocean_slab_mod::snow_wfact
real, parameter snow_wfact
Definition: ocean_slab_mod.F90:59
ocean_slab_mod::tslab
real, dimension(:,:), allocatable, save, public tslab
Definition: ocean_slab_mod.F90:26
dimphy::nslay
integer, save nslay
Definition: dimphy.F90:11
ocean_slab_mod::h_ice_thin
real, parameter h_ice_thin
Definition: ocean_slab_mod.F90:63
ocean_slab_mod::h_ice_max
real, parameter h_ice_max
Definition: ocean_slab_mod.F90:69
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
ocean_slab_mod::sno_den
real, parameter sno_den
Definition: ocean_slab_mod.F90:49
ocean_slab_mod::cyang
real, save, private cyang
Definition: ocean_slab_mod.F90:22
ocean_slab_mod::ocean_slab_ice
subroutine, public ocean_slab_ice(itime, dtime, jour, knon, knindex, tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, AcoefH, AcoefQ, BcoefH, BcoefQ, AcoefU, AcoefV, BcoefU, BcoefV, ps, u1, v1, gustiness, radsol, snow, qsurf, qsol, agesno, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, tsurf_new, dflux_s, dflux_l, swnet)
Definition: ocean_slab_mod.F90:403
indice_sol_mod::epsfra
real, parameter epsfra
Definition: indice_sol_mod.F90:7
ocean_slab_mod::ocean_slab_noice
subroutine, public ocean_slab_noice(itime, dtime, jour, knon, knindex, p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, temp_air, spechum, AcoefH, AcoefQ, BcoefH, BcoefQ, AcoefU, AcoefV, BcoefU, BcoefV, ps, u1, v1, gustiness, tsurf_in, radsol, snow, qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, tsurf_new, dflux_s, dflux_l, qflux)
Definition: ocean_slab_mod.F90:224
ocean_slab_mod::ice_cond
real, parameter ice_cond
Definition: ocean_slab_mod.F90:52
ocean_slab_mod::ice_den
real, parameter ice_den
Definition: ocean_slab_mod.F90:50
ocean_slab_mod::pen_ext
real, parameter pen_ext
Definition: ocean_slab_mod.F90:77
calcul_fluxs_mod
Definition: calcul_fluxs_mod.F90:2
itime
!$Header!integer nvarmx s s itime
Definition: gradsdef.h:20
ocean_slab_mod::snow_min
real, parameter snow_min
Definition: ocean_slab_mod.F90:58
indice_sol_mod::is_sic
integer, parameter is_sic
Definition: indice_sol_mod.F90:3
abort_physic
subroutine abort_physic(modname, message, ierr)
Definition: abort_physic.F90:3
ocean_slab_mod::alb_ice_dry
real, parameter alb_ice_dry
Definition: ocean_slab_mod.F90:74
ocean_slab_mod::h_ice_min
real, parameter h_ice_min
Definition: ocean_slab_mod.F90:62
calcul_fluxs_mod::calcul_fluxs
subroutine calcul_fluxs(knon, nisurf, dtime, tsurf, p1lay, cal, beta, cdragh, cdragq, ps, precip_rain, precip_snow, snow, qsurf, radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, gustiness, fqsat, petAcoef, peqAcoef, petBcoef, peqBcoef, tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
Definition: calcul_fluxs_mod.F90:12
dimphy
Definition: dimphy.F90:1
ocean_slab_mod::sea_den
real, parameter sea_den
Definition: ocean_slab_mod.F90:51
indice_sol_mod::is_oce
integer, parameter is_oce
Definition: indice_sol_mod.F90:3
surface_data
Definition: surface_data.F90:4
indice_sol_mod
Definition: indice_sol_mod.F90:1
ocean_slab_mod::slab_bils
real, dimension(:), allocatable, save, public slab_bils
Definition: ocean_slab_mod.F90:34
ocean_slab_mod::alb_ice_wet
real, parameter alb_ice_wet
Definition: ocean_slab_mod.F90:75