source: LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90 @ 2212

Last change on this file since 2212 was 2209, checked in by Ehouarn Millour, 10 years ago

Update of the slab ocean by Francis Codron. There are now 3 possibilities for the "version_ocean" slab type:
sicOBS = prescribed ice fraction. Water temperature nearby is set to -1.8°C and cannot become lower.
sicNO = ignore sea ice. One can prescribe a fraction, but the nearby ocean evolves freely, depending on surface fluxes: temperature can go below freezing point or above...
sicINT = interactive sea ice.
EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 28.6 KB
RevLine 
[781]1!
2MODULE 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!
[2057]7
8  USE dimphy
9  USE indice_sol_mod
[2209]10  USE surface_data
[2057]11
[781]12  IMPLICIT NONE
[996]13  PRIVATE
[2209]14  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
[781]15
[2209]16!****************************************************************************************
17! Global saved variables
18!****************************************************************************************
19
[2057]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)
[2209]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)
[2057]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)
[2209]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)
[2057]42
[2209]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
[781]81CONTAINS
82!
83!****************************************************************************************
84!
[2057]85  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
86  !, seaice_rst etc
[781]87
[2057]88    use IOIPSL
89
90    INCLUDE "iniprint.h"
91    ! For ok_xxx vars (Ekman...)
92    INCLUDE "clesphys.h"
93
94    ! Input variables
95!****************************************************************************************
96    REAL, INTENT(IN)                         :: dtime
97! Variables read from restart file
98    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
99
100! Local variables
101!****************************************************************************************
102    INTEGER                :: error
103    CHARACTER (len = 80)   :: abort_message
104    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
105
106!****************************************************************************************
107! Allocate surface fraction read from restart file
108!****************************************************************************************
[2209]109    ALLOCATE(fsic(klon), stat = error)
[2057]110    IF (error /= 0) THEN
111       abort_message='Pb allocation tmp_pctsrf_slab'
112       CALL abort_gcm(modname,abort_message,1)
113    ENDIF
[2209]114    fsic(:)=0.
115    WHERE (1.-zmasq(:)>EPSFRA)
116        fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
117    END WHERE
[2057]118
119!****************************************************************************************
[2209]120! Allocate saved variables
[2057]121!****************************************************************************************
[2209]122    ALLOCATE(tslab(klon,nslay), stat=error)
123       IF (error /= 0) CALL abort_gcm &
124         (modname,'pb allocation tslab', 1)
125
[2057]126    ALLOCATE(slab_bils(klon), stat = error)
127    IF (error /= 0) THEN
128       abort_message='Pb allocation slab_bils'
129       CALL abort_gcm(modname,abort_message,1)
130    ENDIF
131    slab_bils(:) = 0.0   
132    ALLOCATE(bils_cum(klon), stat = error)
133    IF (error /= 0) THEN
134       abort_message='Pb allocation slab_bils_cum'
135       CALL abort_gcm(modname,abort_message,1)
136    ENDIF
137    bils_cum(:) = 0.0   
138
[2209]139    IF (version_ocean=='sicINT') THEN
140        ALLOCATE(slab_bilg(klon), stat = error)
141        IF (error /= 0) THEN
142           abort_message='Pb allocation slab_bilg'
143           CALL abort_gcm(modname,abort_message,1)
144        ENDIF
145        slab_bilg(:) = 0.0   
146        ALLOCATE(bilg_cum(klon), stat = error)
147        IF (error /= 0) THEN
148           abort_message='Pb allocation slab_bilg_cum'
149           CALL abort_gcm(modname,abort_message,1)
150        ENDIF
151        bilg_cum(:) = 0.0   
152        ALLOCATE(tice(klon), stat = error)
153        IF (error /= 0) THEN
154           abort_message='Pb allocation slab_tice'
155           CALL abort_gcm(modname,abort_message,1)
156        ENDIF
157        ALLOCATE(seaice(klon), stat = error)
158        IF (error /= 0) THEN
159           abort_message='Pb allocation slab_seaice'
160           CALL abort_gcm(modname,abort_message,1)
161        ENDIF
162    END IF
163
164!****************************************************************************************
165! Define some parameters
166!****************************************************************************************
[2057]167! Layer thickness
168    ALLOCATE(slabh(nslay), stat = error)
169    IF (error /= 0) THEN
170       abort_message='Pb allocation slabh'
171       CALL abort_gcm(modname,abort_message,1)
172    ENDIF
173    slabh(1)=50.
174! cyang = 1/heat capacity of top layer (rho.c.H)
175    cyang=1/(slabh(1)*4.228e+06)
176
177! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
178! pour un calcul à chaque pas de temps, cpl_pas=1
179    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
180    CALL getin('cpl_pas',cpl_pas)
181    print *,'cpl_pas',cpl_pas
[2209]182
[2057]183 END SUBROUTINE ocean_slab_init
184!
185!****************************************************************************************
186!
187  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
188
[996]189    USE limit_read_mod
[1785]190
[996]191!    INCLUDE "clesphys.h"
[781]192
[996]193! Arguments
[781]194!****************************************************************************************
[996]195    INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
196    INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
197    REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
[2057]198    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
[996]199    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
[781]200
201! Local variables
202!****************************************************************************************
203
[2057]204    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
205       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
[996]206    ELSE
[2209]207       pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
208       pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
[2057]209       is_modified=.TRUE.
[996]210    END IF
[781]211
[996]212  END SUBROUTINE ocean_slab_frac
[781]213!
214!****************************************************************************************
215!
216  SUBROUTINE ocean_slab_noice( &
[996]217       itime, dtime, jour, knon, knindex, &
[1067]218       p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
219       AcoefH, AcoefQ, BcoefH, BcoefQ, &
220       AcoefU, AcoefV, BcoefU, BcoefV, &
221       ps, u1, v1, tsurf_in, &
[2209]222       radsol, snow, &
[1067]223       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
[2057]224       tsurf_new, dflux_s, dflux_l, qflux)
[1067]225   
226    USE calcul_fluxs_mod
[1785]227
[781]228    INCLUDE "iniprint.h"
229
230! Input arguments
231!****************************************************************************************
[996]232    INTEGER, INTENT(IN)                  :: itime
233    INTEGER, INTENT(IN)                  :: jour
[781]234    INTEGER, INTENT(IN)                  :: knon
235    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
236    REAL, INTENT(IN)                     :: dtime
237    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
[1067]238    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
[781]239    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
240    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
[1067]241    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
242    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
[781]243    REAL, DIMENSION(klon), INTENT(IN)    :: ps
[1067]244    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1
[996]245    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
[2209]246    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
[781]247
248! In/Output arguments
249!****************************************************************************************
250    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
251   
252! Output arguments
253!****************************************************************************************
254    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
255    REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
[1067]256    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1, flux_v1
[781]257    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
258    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l     
[2057]259    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
[781]260
261! Local variables
262!****************************************************************************************
[2057]263    INTEGER               :: i,ki
[2209]264    ! surface fluxes
[996]265    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
[2209]266    ! flux correction
267    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
268    ! surface wind stress
[1067]269    REAL, DIMENSION(klon) :: u0, v0
270    REAL, DIMENSION(klon) :: u1_lay, v1_lay
[2209]271    ! ice creation
272    REAL                  :: e_freeze, h_new, dfsic
[781]273
274!****************************************************************************************
[996]275! 1) Flux calculation
276!
277!****************************************************************************************
278    cal(:)      = 0.
279    beta(:)     = 1.
280    dif_grnd(:) = 0.
[781]281   
[1067]282! Suppose zero surface speed
283    u0(:)=0.0
284    v0(:)=0.0
285    u1_lay(:) = u1(:) - u0(:)
286    v1_lay(:) = v1(:) - v0(:)
287
[781]288    CALL calcul_fluxs(knon, is_oce, dtime, &
[1067]289         tsurf_in, p1lay, cal, beta, cdragh, ps, &
[781]290         precip_rain, precip_snow, snow, qsurf,  &
291         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
[1067]292         AcoefH, AcoefQ, BcoefH, BcoefQ, &
[781]293         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
294
[1067]295! - Flux calculation at first modele level for U and V
296    CALL calcul_flux_wind(knon, dtime, &
297         u0, v0, u1, v1, cdragm, &
298         AcoefU, AcoefV, BcoefU, BcoefV, &
299         p1lay, temp_air, &
300         flux_u1, flux_v1) 
301
[2057]302! Accumulate total fluxes locally
303    slab_bils(:)=0.
304    DO i=1,knon
305        ki=knindex(i)
[2209]306        slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
307                      -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
[2057]308        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
309! Also taux, tauy, saved vars...
310    END DO
311
[781]312!****************************************************************************************
[2057]313! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
[781]314!
315!****************************************************************************************
[2057]316    lmt_bils(:)=0.
[2209]317    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
318    ! lmt_bils and diff_sst,siv saved by limit_slab
319    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
[2057]320    ! qflux = total QFlux correction (in W/m2)
[781]321
322!****************************************************************************************
[996]323! 3) Recalculate new temperature
[781]324!
[2057]325!***********************************************o*****************************************
326    tsurf_new=tsurf_in
327    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
328        ! Compute transport
329        ! Add read QFlux and SST tendency
330        tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
331        ! Add cumulated surface fluxes
332        tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
333        ! Update surface temperature
334        SELECT CASE(version_ocean)
335        CASE('sicNO')
336            DO i=1,knon
337                ki=knindex(i)
338                tsurf_new(i)=tslab(ki,1)
339            END DO
[2209]340        CASE('sicOBS') ! check for sea ice or tslab below freezing
[2057]341            DO i=1,knon
342                ki=knindex(i)
[2209]343                IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
[2057]344                    tslab(ki,1)=t_freeze
345                END IF
[2209]346                tsurf_new(i)=tslab(ki,1)
[2057]347            END DO
348        CASE('sicINT')
349            DO i=1,knon
350                ki=knindex(i)
[2209]351                IF (fsic(ki).LT.epsfra) THEN ! Free of ice
352                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
353                        ! quantity of new ice formed
354                        e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
355                        ! new ice
356                        tice(ki)=t_freeze
357                        fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
358                        IF (fsic(ki).GT.ice_frac_min) THEN
359                            seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
360                            tslab(ki,1)=t_freeze
361                        ELSE
362                            fsic(ki)=0.
363                        END IF
364                        tsurf_new(i)=t_freeze
365                    ELSE
[2057]366                        tsurf_new(i)=tslab(ki,1)
[2209]367                    END IF
368                ELSE ! ice present
[2057]369                    tsurf_new(i)=t_freeze
[2209]370                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
371                        ! quantity of new ice formed over open ocean
372                        e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
373                                 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
374                        ! new ice height and fraction
375                        h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
376                        dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
377                        h_new=MIN(e_freeze/dfsic,h_ice_max)
378                        ! update tslab to freezing over open ocean only
379                        tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
380                        ! update sea ice
381                        seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
382                                   /(dfsic+fsic(ki))
383                        fsic(ki)=fsic(ki)+dfsic
384                        ! update snow?
385                    END IF !freezing
386                END IF ! sea ice present
[2057]387            END DO
388        END SELECT
389        bils_cum(:)=0.0! clear cumulated fluxes
390    END IF ! coupling time
391  END SUBROUTINE ocean_slab_noice
392!
[781]393!****************************************************************************************
[2209]394
395  SUBROUTINE ocean_slab_ice(   &
396       itime, dtime, jour, knon, knindex, &
397       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
398       AcoefH, AcoefQ, BcoefH, BcoefQ, &
399       AcoefU, AcoefV, BcoefU, BcoefV, &
400       ps, u1, v1, &
401       radsol, snow, qsurf, qsol, agesno, &
402       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
403       tsurf_new, dflux_s, dflux_l, swnet)
404
405   USE calcul_fluxs_mod
406
407   INCLUDE "YOMCST.h"
408
409! Input arguments
[2057]410!****************************************************************************************
[2209]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
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!****************************************************************************************
[2057]458! 1) Flux calculation
459!****************************************************************************************
[2209]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
[2057]498! calcul_fluxs (sens, lat etc)
[2209]499    CALL calcul_fluxs(knon, is_sic, dtime, &
500        tsurf_in, p1lay, cal, beta, cdragh, ps, &
501        precip_rain, precip_snow, snow, qsurf,  &
502        radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
503        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
[2057]509! calcul_flux_wind
[2209]510    CALL calcul_flux_wind(knon, dtime, &
511         u0, v0, u1, v1, cdragm, &
512         AcoefU, AcoefV, BcoefU, BcoefV, &
513         p1lay, temp_air, &
514         flux_u1, flux_v1)
[781]515
[2057]516!****************************************************************************************
[2209]517! 2) Update snow and ice surface
[2057]518!****************************************************************************************
[2209]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
[2057]585! New albedo
[2209]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(:)
[2057]616
617!****************************************************************************************
[2209]618! 3) Recalculate new ocean temperature (add fluxes below ice)
[2057]619!    Melt / freeze from below
620!***********************************************o*****************************************
[2209]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
[2057]681
[2209]682  END SUBROUTINE ocean_slab_ice
[781]683!
684!****************************************************************************************
685!
[2057]686  SUBROUTINE ocean_slab_final
687
688!****************************************************************************************
689! Deallocate module variables
690!****************************************************************************************
691    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
[2209]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)
[2057]698
699  END SUBROUTINE ocean_slab_final
700
[781]701END MODULE ocean_slab_mod
Note: See TracBrowser for help on using the repository browser.