Changeset 2209 for LMDZ5/trunk/libf


Ignore:
Timestamp:
Feb 16, 2015, 8:34:28 AM (10 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

Location:
LMDZ5/trunk/libf/phylmd
Files:
13 edited

Legend:

Unmodified
Added
Removed

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

    r2181 r2209  
    2323
    2424    USE dimphy
    25     USE surface_data, ONLY : type_ocean
     25    USE surface_data, ONLY : type_ocean,version_ocean
    2626    USE limit_read_mod
    2727    USE pbl_surface_mod, ONLY : pbl_surface_newfrac
    2828    USE cpl_mod, ONLY : cpl_receive_frac
    29     USE ocean_slab_mod, ONLY : ocean_slab_frac
     29    USE ocean_slab_mod, ONLY : fsic, ocean_slab_frac
    3030    USE indice_sol_mod
    3131
    … …  
    146146          pctsrf(:,is_sic) = 0.
    147147       END WHERE
     148! Send fractions back to slab ocean if needed
     149       IF (type_ocean == 'slab'.AND. version_ocean.NE.'sicINT') THEN
     150           WHERE (1.-zmasq(:)>EPSFRA)
     151               fsic(:)=pctsrf(:,is_sic)/(1.-zmasq(:))
     152           END WHERE
     153       END IF
    148154
    149155!****************************************************************************************

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

    r1907 r2209  
    9393!****************************************************************************************
    9494
    95     IF (type_ocean == 'couple') THEN
     95IF (type_ocean == 'couple'.OR. &
     96         (type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
    9697       ! limit.nc has not yet been read. Do it now!
    9798       CALL limit_read_tot(itime, dtime, jour, is_modified)

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

    r2057 r2209  
    11! $Header$
    22
    3 SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, diff_sst)
     3SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
    44
    55  USE dimphy
    … …  
    2020  INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
    2121  REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
    22   REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst
     22  REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst, diff_siv
    2323
    2424! Locals variables with attribute SAVE
    2525!****************************************************************************************
    2626  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: bils_save, diff_sst_save
    27 !$OMP THREADPRIVATE(bils_save, diff_sst_save)
     27  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: diff_siv_save
     28!$OMP THREADPRIVATE(bils_save, diff_sst_save, diff_siv_save)
    2829
    2930! Locals variables
    … …  
    3334  INTEGER, DIMENSION(2)    :: start, epais
    3435  REAL, DIMENSION(klon_glo):: bils_glo, sst_l_glo, sst_lp1_glo, diff_sst_glo
     36  REAL, DIMENSION(klon_glo):: siv_l_glo, siv_lp1_glo, diff_siv_glo
    3537  CHARACTER (len = 20)     :: modname = 'limit_slab'
    36   LOGICAL                  :: read_bils,read_sst
     38  LOGICAL                  :: read_bils,read_sst,read_siv
    3739
    3840! End declaration
    … …  
    4951        read_bils=.TRUE.
    5052        read_sst=.TRUE.
     53        read_siv=.TRUE.
    5154       
    5255        ierr = NF90_OPEN ('limit_slab.nc', NF90_NOWRITE, nid)
    … …  
    5457            read_bils=.FALSE.
    5558            read_sst=.FALSE.
     59            read_siv=.FALSE.
    5660        ELSE ! read file
    5761       
    … …  
    6367
    6468!****************************************************************************************
    65 ! 2) Read bils and SST tendency
     69! 2) Read bils and SST/ ice volume tendency
    6670!
    6771!****************************************************************************************
    … …  
    8993        END IF
    9094
     95! Read siv_glo for this day
     96        ierr = NF90_INQ_VARID(nid, 'SICV', nvarid)
     97        IF (ierr /= NF90_NOERR)  THEN
     98            read_siv=.FALSE.
     99        ELSE
     100            start(2) = jour
     101            ierr = NF90_GET_VAR(nid,nvarid,siv_l_glo,start,epais)
     102            IF (ierr /= NF90_NOERR) read_siv=.FALSE.
     103! Read siv_glo for one day ahead
     104            start(2) = jour + 1
     105            IF (start(2) > 360) start(2)=1
     106            ierr = NF90_GET_VAR(nid,nvarid,siv_lp1_glo,start,epais)
     107            IF (ierr /= NF90_NOERR) read_siv=.FALSE.
     108        END IF
     109
    91110!****************************************************************************************
    92111! 5) Close file and distribute variables to all processus
    … …  
    98117        IF (read_sst) THEN
    99118! Calculate difference in temperature between this day and one ahead
    100 !            DO i=1, klon_glo-1
    101 !               diff_sst_glo(i) = sst_lp1_glo(i) - sst_l_glo(i)
    102 !            END DO
    103 !            diff_sst_glo(klon_glo) = sst_lp1_glo(klon_glo) - sst_l_glo(1)
    104119            DO i=1, klon_glo
    105120               diff_sst_glo(i) = sst_lp1_glo(i) - sst_l_glo(i)
    106121            END DO
    107122        END IF !read_sst
     123        IF (read_siv) THEN
     124! Calculate difference in temperature between this day and one ahead
     125            DO i=1, klon_glo
     126               diff_siv_glo(i) = siv_lp1_glo(i) - siv_l_glo(i)
     127            END DO
     128        END IF !read_siv
    108129     ENDIF ! is_mpi_root
    109130
    … …  
    111132       
    112133     IF (.NOT. ALLOCATED(bils_save)) THEN
    113         ALLOCATE(bils_save(klon), diff_sst_save(klon), stat=ierr)
     134        ALLOCATE(bils_save(klon), diff_sst_save(klon), diff_siv_save(klon), stat=ierr)
    114135        IF (ierr /= 0) CALL abort_gcm('limit_slab', 'pb in allocation',1)
    115136     END IF
    116137
    117 ! Giveddefault values if needed
     138! Give default values if needed
    118139     IF (read_bils) THEN
    119140         CALL Scatter(bils_glo, bils_save)
    … …  
    126147         diff_sst_save(:)=0.
    127148     END IF
     149     IF (read_siv) THEN
     150         CALL Scatter(diff_siv_glo, diff_siv_save)
     151     ELSE
     152         diff_siv_save(:)=0.
     153     END IF
    128154     
    129155  ENDIF ! time to read
    … …  
    131157  lmt_bils(:) = bils_save(:)
    132158  diff_sst(:) = diff_sst_save(:)
     159  diff_siv(:) = diff_siv_save(:)
     160
    133161 
    134162END SUBROUTINE limit_slab

  • 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

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

    r2189 r2209  
    3333  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
    3434  !$OMP THREADPRIVATE(fder)
    35   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: snow   ! snow at surface
     35  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: snow   ! snow at surface
    3636  !$OMP THREADPRIVATE(snow)
    3737  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface

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

    r2198 r2209  
    238238    rugos = rugos_omp
    239239    WRITE (*, *) 'iniaqua: rugos=', rugos
    240     zmasq(:) = pctsrf(:, is_oce)
     240    zmasq(:) = pctsrf(:, is_ter)
    241241
    242242    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)

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

    r2190 r2209  
    88  USE fonte_neige_mod,  ONLY : fonte_neige_init
    99  USE pbl_surface_mod,  ONLY : pbl_surface_init
    10   USE surface_data,     ONLY : type_ocean
     10  USE surface_data,     ONLY : type_ocean, version_ocean
    1111  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
    1212       du_gwd_rando, dv_gwd_rando, entr_therm, f0, falb1, falb2, fm_therm, &
    … …  
    2323  USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
    2424  USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
    25   USE ocean_slab_mod, ONLY: tslab, ocean_slab_init
    26 
     25  USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init
    2726
    2827  IMPLICIT none
    … …  
    11611160  ! Initialize Slab variables
    11621161  IF ( type_ocean == 'slab' ) THEN
    1163       ALLOCATE(tslab(klon,nslay), stat=ierr)
    1164       IF (ierr /= 0) CALL abort_gcm &
    1165          ('phyetat0', 'pb allocation tslab', 1)
     1162      print*, "calling slab_init"
     1163      CALL ocean_slab_init(dtime, pctsrf)
     1164      ! tslab
    11661165      CALL get_field("tslab", tslab, found)
    11671166      IF (.NOT. found) THEN
    … …  
    11691168          PRINT*, "Initialisation a tsol_oce"
    11701169          DO i=1,nslay
    1171               tslab(:,i)=ftsol(:,is_oce)
     1170              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
    11721171          END DO
    11731172      END IF
    1174       print*, "calling slab_init"
    1175       CALL ocean_slab_init(dtime, pctsrf)
     1173      ! Sea ice variables
     1174      IF (version_ocean == 'sicINT') THEN
     1175          CALL get_field("slab_tice", tice, found)
     1176          IF (.NOT. found) THEN
     1177              PRINT*, "phyetat0: Le champ <tice> est absent"
     1178              PRINT*, "Initialisation a tsol_sic"
     1179                  tice(:)=ftsol(:,is_sic)
     1180          END IF
     1181          CALL get_field("seaice", seaice, found)
     1182          IF (.NOT. found) THEN
     1183              PRINT*, "phyetat0: Le champ <seaice> est absent"
     1184              PRINT*, "Initialisation a 0/1m suivant fraction glace"
     1185              seaice(:)=0.
     1186              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
     1187                  seaice=917.
     1188              END WHERE
     1189          END IF
     1190      END IF !sea ice INT
    11761191  END IF ! Slab       
    11771192

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

    r2188 r2209  
    1616  USE indice_sol_mod
    1717  USE surface_data
    18   USE ocean_slab_mod, ONLY : tslab
     18  USE ocean_slab_mod, ONLY : tslab, seaice, tice, fsic
    1919
    2020  IMPLICIT none
    … …  
    106106
    107107  ! BP ajout des fraction de chaque sous-surface
     108
     109  ! Get last fractions from slab ocean
     110  IF (type_ocean == 'slab' .AND. version_ocean == "sicINT") THEN
     111      WHERE (1.-zmasq(:).GT.EPSFRA)
     112          pctsrf(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
     113          pctsrf(:,is_sic)=fsic(:)*(1.-zmasq(:))
     114      END WHERE
     115  END IF
    108116
    109117  ! 1. fraction de terre
    … …  
    352360  IF (type_ocean == 'slab') THEN
    353361      CALL put_field("tslab", "Slab ocean temperature", tslab)
     362      IF (version_ocean == 'sicINT') THEN
     363          CALL put_field("seaice", "Slab seaice (kg/m2)", seaice)
     364          CALL put_field("slab_tice", "Slab sea ice temperature", tice)
     365      END IF
    354366  END IF
    355367

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

    r2183 r2209  
    465465  TYPE(ctrl_out), SAVE :: o_slab_bils = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
    466466    'slab_bils', 'flux atmos - slab ponderes foce', 'W/m2', (/ ('', i=1, 9) /))
     467  TYPE(ctrl_out), SAVE :: o_slab_bilg = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     468    'slab_bilg', 'flux glace - slab ponderes fsic', 'W/m2', (/ ('', i=1, 9) /))
    467469  TYPE(ctrl_out), SAVE :: o_slab_qflux = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
    468470    'slab_qflux', 'Correction flux slab', 'W/m2', (/ ('', i=1, 9) /))
    469471  TYPE(ctrl_out), SAVE :: o_tslab = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
    470472    'tslab', 'Temperature ocean slab', 'K', (/ ('', i=1, 9) /))
     473  TYPE(ctrl_out), SAVE :: o_slab_tice = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     474    'slab_tice', 'Temperature banquise slab', 'K', (/ ('', i=1, 9) /))
     475  TYPE(ctrl_out), SAVE :: o_slab_sic = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     476    'seaice', 'Epaisseur banquise slab', 'kg/m2', (/ ('', i=1, 9) /))
    471477  TYPE(ctrl_out), SAVE :: o_ale_bl = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11 /), &
    472478    'ale_bl', 'ALE BL', 'm2/s2', (/ ('', i=1, 9) /))
    … …  
    10081014      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'ages_oce',"Snow age", "day", (/ ('', i=1, 9) /)), &
    10091015      ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'ages_sic',"Snow age", "day", (/ ('', i=1, 9) /)) /)
     1016
     1017  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_snow_srf     = (/ &
     1018      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_ter', "Snow", "kg/m2", (/ ('', i=1, 9) /)), &
     1019      ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_lic', "Snow", "kg/m2", (/ ('', i=1, 9) /)), &
     1020      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_oce',"Snow", "kg/m2", (/ ('', i=1, 9) /)), &
     1021      ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_sic',"Snow", "kg/m2", (/ ('', i=1, 9) /)) /)
    10101022
    10111023  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_rugs_srf     = (/ &

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

    r2194 r2209  
    8080         o_alp_bl_conv, o_alp_bl_stat, &
    8181         o_slab_qflux, o_tslab, o_slab_bils, &
     82         o_slab_bilg, o_slab_sic, o_slab_tice, &
    8283         o_weakinv, o_dthmin, o_cldtau, &
    8384         o_cldemi, o_pr_con_l, o_pr_con_i, &
    … …  
    113114         o_rnebls, o_rhum, o_ozone, o_ozone_light, &
    114115         o_dtphy, o_dqphy, o_albe_srf, o_rugs_srf, &
    115          o_ages_srf, o_alb1, o_alb2, o_tke, &
     116         o_ages_srf, o_snow_srf, o_alb1, o_alb2, o_tke, &
    116117         o_tke_max, o_kz, o_kz_max, o_clwcon, &
    117118         o_dtdyn, o_dqdyn, o_dudyn, o_dvdyn, &
    … …  
    227228         bils_ec,bils_ech, bils_tke, bils_kinetic, bils_latent, bils_enthalp, &
    228229         itau_con, nfiles, clef_files, nid_files, zvstr_gwd_rando
    229     USE ocean_slab_mod, only: tslab, slab_bils
     230    USE ocean_slab_mod, only: tslab, slab_bils, slab_bilg, tice, seaice
     231    USE pbl_surface_mod, only: snow
    230232    USE indice_sol_mod, only: nbsrf
    231233    USE infotrac, only: nqtot, nqo, type_trac
    232234    USE comgeomphy, only: airephy
    233     USE surface_data, only: type_ocean, ok_veget, ok_snow
     235    USE surface_data, only: type_ocean, version_ocean, ok_veget, ok_snow
    234236!    USE aero_mod, only: naero_spc
    235237    USE aero_mod, only: naero_tot, id_STRAT_phy
    … …  
    400402       CALL histwrite_phy(o_pluc, zx_tmp_fi2d)
    401403       CALL histwrite_phy(o_snow, snow_fall)
    402        CALL histwrite_phy(o_msnow, snow_o)
     404       CALL histwrite_phy(o_msnow, zxsnow)
    403405       CALL histwrite_phy(o_fsnow, zfra_o)
    404406       CALL histwrite_phy(o_evap, evap)
    … …  
    514516
    515517       IF (ok_snow) THEN
    516           CALL histwrite_phy(o_snowsrf, zxsnow)
     518          CALL histwrite_phy(o_snowsrf, snow_o)
    517519          CALL histwrite_phy(o_qsnow, qsnow)
    518520          CALL histwrite_phy(o_snowhgt,snowhgt)
    … …  
    755757              CALL histwrite_phy(o_tslab, tslab)
    756758          END IF
     759          IF (version_ocean=='sicINT') THEN
     760              CALL histwrite_phy(o_slab_bilg, slab_bilg)
     761              CALL histwrite_phy(o_slab_tice, tice)
     762              CALL histwrite_phy(o_slab_sic, seaice)
     763          END IF
    757764       ENDIF !type_ocean == force/slab
    758765       CALL histwrite_phy(o_weakinv, weak_inversion)
    … …  
    970977          IF (vars_defined) zx_tmp_fi2d(1 : klon) = agesno( 1 : klon, nsrf)
    971978          CALL histwrite_phy(o_ages_srf(nsrf), zx_tmp_fi2d)
     979          IF (vars_defined) zx_tmp_fi2d(1 : klon) = snow( 1 : klon, nsrf)
     980          CALL histwrite_phy(o_snow_srf(nsrf), zx_tmp_fi2d)
    972981       ENDDO !nsrf=1, nbsrf
    973982       CALL histwrite_phy(o_alb1, albsol1)

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

    r2178 r2209  
    123123            AcoefU, AcoefV, BcoefU, BcoefV, &
    124124            ps, u1, v1, tsurf_in, &
    125             radsol, snow, agesno, &
     125            radsol, snow, &
    126126            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    127127            tsurf_new, dflux_s, dflux_l, lmt_bils)

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

    r2057 r2209  
    2525  USE ocean_forced_mod, ONLY : ocean_forced_ice
    2626  USE ocean_cpl_mod, ONLY    : ocean_cpl_ice
     27  USE ocean_slab_mod, ONLY   : ocean_slab_ice
    2728  USE indice_sol_mod
    2829
    … …  
    108109            tsurf_new, dflux_s, dflux_l)
    109110       
    110 !    ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN
    111 !       CALL ocean_slab_ice( &
    112 !          itime, dtime, jour, knon, knindex, &
    113 !          debut, &
    114 !          tsurf, p1lay, cdragh, precip_rain, precip_snow, temp_air, spechum,&
    115 !          AcoefH, AcoefQ, BcoefH, BcoefQ, &
    116 !          ps, u1, v1, &
    117 !          radsol, snow, qsurf, qsol, agesno, tsoil, &
    118 !          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
    119 !          tsurf_new, dflux_s, dflux_l)
    120 !
     111    ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN
     112       CALL ocean_slab_ice( &
     113          itime, dtime, jour, knon, knindex, &
     114          tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum,&
     115          AcoefH, AcoefQ, BcoefH, BcoefQ, &
     116            AcoefU, AcoefV, BcoefU, BcoefV, &
     117          ps, u1, v1, &
     118          radsol, snow, qsurf, qsol, agesno, &
     119          alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
     120          tsurf_new, dflux_s, dflux_l, swnet)
     121
    121122      ELSE ! type_ocean=force or slab +sicOBS or sicNO
    122123       CALL ocean_forced_ice( &

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

    r2075 r2209  
    77  REAL, PARAMETER        :: tau_gl=86400.*5.
    88  REAL, PARAMETER        :: calsno=1./(2.3867e+06*.15)
    9   REAL, PARAMETER        :: t_freeze=271.35
    10   REAL, PARAMETER        :: t_melt=273.15
    119 
    1210  LOGICAL, SAVE          :: ok_veget      ! true for use of vegetation model ORCHIDEE
Note: See TracChangeset for help on using the changeset viewer.