Ignore:
Timestamp:
Feb 16, 2015, 8:34:28 AM (9 years ago)
Author:
Ehouarn Millour
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90

    r2057 r2209  
    88  USE dimphy
    99  USE indice_sol_mod
     10  USE surface_data
    1011
    1112  IMPLICIT NONE
    1213  PRIVATE
    13   PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice!, ocean_slab_ice
     14  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
     15
     16!****************************************************************************************
     17! Global saved variables
     18!****************************************************************************************
    1419
    1520  INTEGER, PRIVATE, SAVE                           :: cpl_pas
    … …  
    2126  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
    2227  !$OMP THREADPRIVATE(tslab)
    23   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE  :: pctsrf
    24   !$OMP THREADPRIVATE(pctsrf)
     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)
    2534  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
    2635  !$OMP THREADPRIVATE(slab_bils)
    2736  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
    2837  !$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!****************************************************************************************
    2980
    3081CONTAINS
    … …  
    56107! Allocate surface fraction read from restart file
    57108!****************************************************************************************
    58     ALLOCATE(pctsrf(klon,nbsrf), stat = error)
     109    ALLOCATE(fsic(klon), stat = error)
    59110    IF (error /= 0) THEN
    60111       abort_message='Pb allocation tmp_pctsrf_slab'
    61112       CALL abort_gcm(modname,abort_message,1)
    62113    ENDIF
    63     pctsrf(:,:) = pctsrf_rst(:,:)
    64 
    65 !****************************************************************************************
    66 ! Allocate local variables
    67 !****************************************************************************************
     114    fsic(:)=0.
     115    WHERE (1.-zmasq(:)>EPSFRA)
     116        fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
     117    END WHERE
     118
     119!****************************************************************************************
     120! Allocate saved variables
     121!****************************************************************************************
     122    ALLOCATE(tslab(klon,nslay), stat=error)
     123       IF (error /= 0) CALL abort_gcm &
     124         (modname,'pb allocation tslab', 1)
     125
    68126    ALLOCATE(slab_bils(klon), stat = error)
    69127    IF (error /= 0) THEN
    … …  
    79137    bils_cum(:) = 0.0   
    80138
     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!****************************************************************************************
    81167! Layer thickness
    82168    ALLOCATE(slabh(nslay), stat = error)
    … …  
    94180    CALL getin('cpl_pas',cpl_pas)
    95181    print *,'cpl_pas',cpl_pas
     182
    96183 END SUBROUTINE ocean_slab_init
    97184!
    … …  
    101188
    102189    USE limit_read_mod
    103     USE surface_data
    104190
    105191!    INCLUDE "clesphys.h"
    … …  
    119205       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
    120206    ELSE
    121        pctsrf_chg(:,:)=pctsrf(:,:)
     207       pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
     208       pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
    122209       is_modified=.TRUE.
    123210    END IF
    … …  
    133220       AcoefU, AcoefV, BcoefU, BcoefV, &
    134221       ps, u1, v1, tsurf_in, &
    135        radsol, snow, agesno, &
     222       radsol, snow, &
    136223       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    137224       tsurf_new, dflux_s, dflux_l, qflux)
    138225   
    139226    USE calcul_fluxs_mod
    140     USE surface_data
    141227
    142228    INCLUDE "iniprint.h"
    … …  
    158244    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1
    159245    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
     246    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
    160247
    161248! In/Output arguments
    162249!****************************************************************************************
    163     REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
    164250    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
    165     REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
    166251   
    167252! Output arguments
    … …  
    177262!****************************************************************************************
    178263    INTEGER               :: i,ki
     264    ! surface fluxes
    179265    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
    180     REAL, DIMENSION(klon) :: diff_sst, lmt_bils
     266    ! flux correction
     267    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
     268    ! surface wind stress
    181269    REAL, DIMENSION(klon) :: u0, v0
    182270    REAL, DIMENSION(klon) :: u1_lay, v1_lay
     271    ! ice creation
     272    REAL                  :: e_freeze, h_new, dfsic
    183273
    184274!****************************************************************************************
    … …  
    189279    beta(:)     = 1.
    190280    dif_grnd(:) = 0.
    191     agesno(:)   = 0.
    192281   
    193282! Suppose zero surface speed
    … …  
    215304    DO i=1,knon
    216305        ki=knindex(i)
    217         slab_bils(ki)=(fluxlat(i)+fluxsens(i)+radsol(i))*pctsrf(ki,is_oce)/(1.-zmasq(ki))
     306        slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
     307                      -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
    218308        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
    219309! Also taux, tauy, saved vars...
    … …  
    225315!****************************************************************************************
    226316    lmt_bils(:)=0.
    227     CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst) ! global pour un processus
    228     ! lmt_bils and diff_sst saved by limit_slab
    229     qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.
     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.
    230320    ! qflux = total QFlux correction (in W/m2)
    231321
    … …  
    248338                tsurf_new(i)=tslab(ki,1)
    249339            END DO
    250         CASE('sicOBS') ! check for sea ice or tsurf below freezing
     340        CASE('sicOBS') ! check for sea ice or tslab below freezing
    251341            DO i=1,knon
    252342                ki=knindex(i)
    253                 IF ((tslab(ki,1).LT.t_freeze).OR.(pctsrf(ki,is_sic).GT.epsfra)) THEN
    254                     tsurf_new(i)=t_freeze
     343                IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
    255344                    tslab(ki,1)=t_freeze
    256                 ELSE
    257                     tsurf_new(i)=tslab(ki,1)
    258345                END IF
     346                tsurf_new(i)=tslab(ki,1)
    259347            END DO
    260348        CASE('sicINT')
    261349            DO i=1,knon
    262350                ki=knindex(i)
    263                 IF (pctsrf(ki,is_sic).LT.epsfra) THEN ! Free of ice
    264                     IF (tslab(ki,1).GT.t_freeze) THEN
     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
    265366                        tsurf_new(i)=tslab(ki,1)
    266                     ELSE
    267                         tsurf_new(i)=t_freeze
    268                         ! Call new ice routine
    269                         tslab(ki,1)=t_freeze
    270                     END IF
    271                 ELSE ! ice present, tslab update completed in slab_ice
     367                    END IF
     368                ELSE ! ice present
    272369                    tsurf_new(i)=t_freeze
    273                 END IF !ice free
     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
    274387            END DO
    275388        END SELECT
    … …  
    279392!
    280393!****************************************************************************************
    281 !
    282 !  SUBROUTINE ocean_slab_ice(   &
    283 !       itime, dtime, jour, knon, knindex, &
    284 !       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
    285 !       AcoefH, AcoefQ, BcoefH, BcoefQ, &
    286 !       AcoefU, AcoefV, BcoefU, BcoefV, &
    287 !       ps, u1, v1, &
    288 !       radsol, snow, qsurf, qsol, agesno, tsoil, &
    289 !       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    290 !       tsurf_new, dflux_s, dflux_l)
    291 !
     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
     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
     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
    292457!****************************************************************************************
    293458! 1) Flux calculation
    294459!****************************************************************************************
    295 ! set beta, cal etc. depends snow / ice surf ?
     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
    296498! calcul_fluxs (sens, lat etc)
     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
    297509! calcul_flux_wind
    298 
    299 !****************************************************************************************
    300 ! 2) Update surface
    301 !****************************************************************************************
    302 ! neige, fonte
    303 ! flux glace-ocean
    304 ! update temperature
    305 ! neige precip, evap
    306 ! Melt snow & ice from above
     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)
     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
    307585! New albedo
    308 
    309 !****************************************************************************************
    310 ! 3) Recalculate new ocean temperature
     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)
    311619!    Melt / freeze from below
    312620!***********************************************o*****************************************
    313 
    314 
    315 !  END SUBROUTINE ocean_slab_ice
     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
    316683!
    317684!****************************************************************************************
    318685!
    319686  SUBROUTINE ocean_slab_final
    320   !, seaice_rst etc
    321 
    322     ! For ok_xxx vars (Ekman...)
    323     INCLUDE "clesphys.h"
    324687
    325688!****************************************************************************************
    326689! Deallocate module variables
    327 !
    328 !****************************************************************************************
    329     IF (ALLOCATED(pctsrf)) DEALLOCATE(pctsrf)
     690!****************************************************************************************
    330691    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)
    331698
    332699  END SUBROUTINE ocean_slab_final
Note: See TracChangeset for help on using the changeset viewer.