Changeset 5662 for LMDZ6/trunk/libf


Ignore:
Timestamp:
May 20, 2025, 4:24:41 PM (3 weeks ago)
Author:
Laurent Fairhead
Message:

Ajout du modèle thermodynamique de glace de mer interactive améliorant les flux échangés à la surface de la banquise (Doctorat de Nicolas Michalezyk, Contact : Nicolas Michaleyk, Guillaume Gastineau)

Location:
LMDZ6/trunk/libf
Files:
24 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90

    r5296 r5662  
    110110  REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp,phiso
    111111  REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
     112  !GG
     113  REAL, DIMENSION(klon)               :: hice, tsic, bilg_cumul
     114  !GG
    112115  REAL, DIMENSION(klon,nbsrf)         :: qsurf, snsrf
    113116  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
    … …  
    142145                   iflag_cldcon,                                        &
    143146                   ratqsbas,ratqshaut,tau_ratqs,            &
     147             !GG      iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
    144148                   ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat,     &
    145149                   aerosol_couple, chemistry_couple, flag_aerosol, flag_aerosol_strat,  &
    … …  
    278282  detr_therm       = 0.
    279283  awake_s = 0.
     284  !GG
     285  hice             = 1.0
     286  tsic             = tsol
     287  bilg_cumul       = 0.
     288  !GG
    280289
    281290  CALL fonte_neige_init(run_off_lic_0)
    282   CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )
     291!GG  CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )
     292  CALL pbl_surface_init( fder, snsrf, qsurf, tsoil, hice, tsic, bilg_cumul)
     293  !GG
    283294
    284295  IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) then

  • LMDZ6/trunk/libf/phylmd/Dust/phys_output_write_spl_mod.F90

    r5337 r5662  
    378378         scdnc, cldncl, reffclws, reffclwc, cldnvi, &
    379379         lcc, lcc3d, lcc3dcon, lcc3dstra, reffclwtop
    380     USE ocean_slab_mod, ONLY: tslab, slab_bilg, tice, seaice
     380    USE ocean_slab_mod, ONLY: tslab, slab_bilg, tice_slab, seaice
    381381    USE pbl_surface_mod, ONLY: snow
    382382    USE indice_sol_mod, ONLY: nbsrf
    … …  
    955955          IF (version_ocean=='sicINT') THEN
    956956              CALL histwrite_phy(o_slab_bilg, slab_bilg)
    957               CALL histwrite_phy(o_slab_tice, tice)
     957              CALL histwrite_phy(o_slab_tice, tice_slab)
    958958              CALL histwrite_phy(o_slab_sic, seaice)
    959959          ENDIF

  • LMDZ6/trunk/libf/phylmd/conf_phys_m.f90

    r5651 r5662  
    219219    INTEGER, SAVE :: iflag_cycle_diurne_omp
    220220    LOGICAL, SAVE :: soil_model_omp,liqice_in_radocond_omp
     221    !GG
     222    INTEGER,SAVE :: iflag_seaice_omp, iflag_seaice_alb_omp
     223    INTEGER,SAVE :: iflag_leads_omp
     224    REAL,SAVE :: sice_cond_omp, sisno_cond_omp
     225    REAL,SAVE :: sisno_den_omp, sisno_min_omp, sithick_min_omp
     226    REAL,SAVE :: sisno_wfact_omp, amax_n_omp, amax_s_omp, rn_alb_sdry_omp
     227    REAL,SAVE :: rn_alb_smlt_omp, rn_alb_idry_omp, rn_alb_imlt_omp
     228    REAL,SAVE :: si_pen_frac_omp,si_pen_ext_omp
     229    REAL,SAVE :: fseaN_omp, fseaS_omp
     230    !GG
    221231    LOGICAL, SAVE :: ok_orodr_omp, ok_orolf_omp, ok_limitvrai_omp
    222232    INTEGER, SAVE :: nbapp_rad_omp, iflag_con_omp
    … …  
    879889    CALL getin('liqice_in_radocond',liqice_in_radocond_omp)
    880890
     891    !GG
     892    !Config  Key  = iflag_seaice
     893    !Config  Desc = Flag de sea ice modele
     894    !Config  Def  = 0
     895    !Config  Help = Flag  pour sea ice thermo  les options suivantes existent :
     896    !Config         0 pour defaut,
     897    !Config         1 pour lu avec limit.nc,
     898    !Config         2 pour interactive sic
     899    iflag_seaice_omp = 0
     900    CALL getin('iflag_seaice',iflag_seaice_omp)
     901
     902    !Config  Key  = iflag_seaice_alb
     903    !Config  Desc = Flag de sea ice thickness
     904    !Config  Def  = 0
     905    !Config  Help = Flag  pour albedo sea ice les options suivantes existent :
     906    !Config         0 pour defaut,
     907    !Config         1 pour new,
     908    !Config         2 pour LIM3,
     909    iflag_seaice_alb_omp = 0
     910    CALL getin('iflag_seaice_alb',iflag_seaice_alb_omp)
     911
     912    !Config  Key  = iflag_seaice_alb
     913    !Config  Desc = Flag de sea ice leads
     914    !Config  Def  = 0
     915    !Config  Help = Flag  pour les leads les options suivantes existent :
     916    !Config         0 pour defaut,
     917    !Config         1 pour correction sur leads,
     918    iflag_leads_omp = 0
     919    CALL getin('iflag_leads',iflag_leads_omp)
     920
     921    !Config key  = sice_cond_omp
     922    !Config Desc = conductivity of ice
     923    !Config Def  = 2.17
     924    !Config Help = pour le modele de glace de mer
     925    sice_cond_omp = 2.17
     926    CALL getin('sice_cond', sice_cond_omp)
     927
     928    !Config key  = sisno_cond_omp
     929    !Config Desc = conductivity of snow
     930    !Config Def  = 0.31
     931    !Config Help = pour le modele de glace de mer
     932    sisno_cond_omp = 0.31
     933    CALL getin('sno_cond', sisno_cond_omp)
     934
     935    !Config key  = sisno_den_omp
     936    !Config Desc = density of snow
     937    !Config Def  = 300
     938    !Config Help = pour le modele de glace de mer
     939    sisno_den_omp = 300
     940    CALL getin('sisno_den', sisno_den_omp)
     941
     942    !Config key  = sisno_min
     943    !Config Desc = minimum snow thickness over sea ice
     944    !Config Def  = 0.05
     945    !Config Help = pour le modele de glace de mer
     946    sisno_min_omp = 0.05
     947    CALL getin('sisno_min', sisno_min_omp)
     948
     949    !Config key  = sithick_min
     950    !Config Desc = minimum sea ice thickness
     951    !Config Def  = 0.10
     952    !Config Help = pour le modele de glace de mer
     953    sithick_min_omp = 0.10
     954    CALL getin('sithick_min', sithick_min_omp)
     955
     956    !Config key  = sisno_wfact
     957    !Config Desc = max fraction of falling snow blown into ocean
     958    !Config Def  = 0.4
     959    !Config Help = pour le modele de glace de mer
     960    sisno_wfact_omp = 0.4
     961    CALL getin('sisno_wfact', sisno_wfact_omp)
     962
     963    !Config key  = amax_n
     964    !Config Desc = leads fraction in NH given as a maximum sea ice concentration
     965    !Config Def  = 0.997
     966    !Config Help = pour le modele de glace de mer
     967    amax_n_omp = 0.997
     968    CALL getin('amax_n', amax_n_omp)
     969
     970    !Config key  = amax_s
     971    !Config Desc = leads fraction in SH given as a maximum sea ice concentration
     972    !Config Def  = 0.95
     973    !Config Help = pour le modele de glace de mer
     974    amax_s_omp = 0.95
     975    CALL getin('amax_s', amax_s_omp)
     976
     977    !Config key  = rn_alb_sdry
     978    !Config Desc = dry snow over sea ice albedo
     979    !Config Def  = 0.85
     980    !Config Help = pour le modele de glace de mer
     981    rn_alb_sdry_omp = 0.85
     982    CALL getin('rn_alb_sdry', rn_alb_sdry_omp)
     983
     984    !Config key  = rn_alb_smlt
     985    !Config Desc = wet snow over sea ice albedo
     986    !Config Def  = 0.55
     987    !Config Help = pour le modele de glace de mer
     988    rn_alb_smlt_omp = 0.55
     989    CALL getin('rn_alb_smlt', rn_alb_smlt_omp)
     990
     991    !Config key  = rn_alb_idry
     992    !Config Desc = dry sea ice albedo
     993    !Config Def  = 0.75
     994    !Config Help = pour le modele de glace de mer
     995    rn_alb_idry_omp = 0.75
     996    CALL getin('rn_alb_idry', rn_alb_idry_omp)
     997
     998    !Config key  = rn_alb_imlt
     999    !Config Desc = wet sea ice albedo
     1000    !Config Def  = 0.66
     1001    !Config Help = pour le modele de glace de mer
     1002    rn_alb_imlt_omp = 0.66
     1003    CALL getin('rn_alb_imlt', rn_alb_imlt_omp)
     1004
     1005    !Config key  = si_pen_frac
     1006    !Config Desc = fraction of shortwave penetrating into the ice
     1007    !Config Def  = 0.3
     1008    !Config Help = pour le modele de glace de mer
     1009    si_pen_frac_omp = 0.3
     1010    CALL getin('si_pen_frac', si_pen_frac_omp)
     1011
     1012    !Config key  = si_pen_ext
     1013    !Config Desc = extinction length of penetrating shortwave (m-1)
     1014    !Config Def  = 1.5
     1015    !Config Help = pour le modele de glace de mer
     1016    si_pen_ext_omp =1.5
     1017    CALL getin('si_pen_ext', si_pen_ext_omp)
     1018
     1019    !Config key  = fseaN
     1020    !Config Desc = HF from ocean below ice Northern Hemisphere
     1021    !Config Def  = 2.0
     1022    !Config Help = pour le modele de glace de mer
     1023    fseaN_omp =2.0
     1024    CALL getin('fseaN', fseaN_omp)
     1025
     1026    !Config key  = fseaS
     1027    !Config Desc = HF from ocean below ice Southern Hemisphere
     1028    !Config Def  = 4.0
     1029    !Config Help = pour le modele de glace de mer
     1030    fseaS_omp =4.0
     1031    CALL getin('fseaS', fseaS_omp)
     1032
     1033    !GG
     1034
    8811035    !Config  Key  = ok_orodr
    8821036    !Config  Desc = Orodr ???
    … …  
    23442498    soil_model = soil_model_omp
    23452499    liqice_in_radocond = liqice_in_radocond_omp
     2500! GG
     2501    sice_cond = sice_cond_omp
     2502    sisno_cond = sisno_cond_omp
     2503    iflag_seaice = iflag_seaice_omp
     2504    iflag_seaice_alb = iflag_seaice_alb_omp
     2505    iflag_leads = iflag_leads_omp
     2506    sisno_den = sisno_den_omp
     2507    sisno_min = sisno_min_omp
     2508    sithick_min = sithick_min_omp
     2509    sisno_wfact = sisno_wfact_omp
     2510    amax_n = amax_n_omp
     2511    amax_s = amax_s_omp
     2512    rn_alb_sdry = rn_alb_sdry_omp
     2513    rn_alb_smlt = rn_alb_smlt_omp
     2514    rn_alb_idry = rn_alb_idry_omp
     2515    rn_alb_imlt = rn_alb_imlt_omp
     2516    fseaN = fseaN_omp
     2517    fseaS = fseaS_omp
     2518    si_pen_ext = si_pen_ext_omp
     2519    si_pen_frac = si_pen_frac_omp
     2520! GG
     2521
    23462522    ok_orodr = ok_orodr_omp
    23472523    ok_orolf = ok_orolf_omp
    … …  
    29363112    WRITE(lunout,*) ' var_fco2_land_cor = ', var_fco2_land_cor
    29373113    WRITE(lunout,*) ' dms_cycle_cpl = ', dms_cycle_cpl
     3114    !GG
     3115    WRITE(lunout,*) ' iflag_seaice = ', iflag_seaice
     3116    WRITE(lunout,*) ' iflag_seaice_alb = ', iflag_seaice_alb
     3117    WRITE(lunout,*) ' iflag_leads = ', iflag_leads
     3118    WRITE(lunout,*) ' sice_cond = ', sice_cond
     3119    WRITE(lunout,*) ' sisno_cond = ', sisno_cond
     3120    WRITE(lunout,*) ' sisno_den = ', sisno_den
     3121    WRITE(lunout,*) ' sisno_min = ', sisno_min
     3122    WRITE(lunout,*) ' sithick_min = ', sithick_min
     3123    WRITE(lunout,*) ' sisno_wfact = ', sisno_wfact
     3124    WRITE(lunout,*) ' amax_n = ', amax_n
     3125    WRITE(lunout,*) ' amax_s = ', amax_s
     3126    WRITE(lunout,*) ' rn_alb_sdry = ', rn_alb_sdry
     3127    WRITE(lunout,*) ' rn_alb_smlt = ', rn_alb_smlt
     3128    WRITE(lunout,*) ' rn_alb_idry = ', rn_alb_idry
     3129    WRITE(lunout,*) ' rn_alb_imlt = ', rn_alb_imlt
     3130    WRITE(lunout,*) ' si_pen_frac = ', si_pen_frac
     3131    WRITE(lunout,*) ' si_pen_ext = ', si_pen_ext
     3132    WRITE(lunout,*) ' fseaN = ', fseaN
     3133    WRITE(lunout,*) ' fseaS = ', fseaS
     3134!GG
    29383135    WRITE(lunout,*) ' n2o_cycle_cpl = ', n2o_cycle_cpl
    29393136    WRITE(lunout,*) ' iflag_tsurf_inlandsis = ', iflag_tsurf_inlandsis

  • LMDZ6/trunk/libf/phylmd/create_etat0_unstruct_mod.f90

    r5296 r5662  
    305305
    306306    CALL fonte_neige_init(run_off_lic_0)
    307     CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )
     307    !GG
     308    ! CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil )
     309    CALL pbl_surface_init( fder, snsrf, qsurf, tsoil, hice, tice, bilg_cumul )
     310    !GG
     311
    308312
    309313    IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) then

  • LMDZ6/trunk/libf/phylmd/limit_read_mod.f90

    r5268 r5662  
    2121!$OMP THREADPRIVATE(albedo) 
    2222  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: sst
    23 !$OMP THREADPRIVATE(sst) 
     23!$OMP THREADPRIVATE(sst)
     24!GG
     25  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: sih
     26!$OMP THREADPRIVATE(sih)
     27!GG
    2428  LOGICAL,SAVE :: read_continents=.FALSE.
    2529!$OMP THREADPRIVATE(read_continents)
    … …  
    145149  END SUBROUTINE limit_read_sst
    146150
     151!GG
     152  SUBROUTINE limit_read_hice(knon, knindex, hice_out)
     153!
     154! This subroutine returns the sea surface temperature already read from limit.nc.
     155!
     156    USE dimphy, ONLY : klon
     157
     158    INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
     159    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
     160    REAL, DIMENSION(klon), INTENT(OUT)   :: hice_out
     161
     162    INTEGER :: i
     163
     164    DO i = 1, knon
     165       hice_out(i) = sih(knindex(i))
     166    END DO
     167
     168  END SUBROUTINE limit_read_hice
     169!GG
    147170!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    148171!!
    … …  
    164187    USE mod_grid_phy_lmdz
    165188    USE mod_phys_lmdz_para
    166     USE surface_data, ONLY : type_ocean, ok_veget
     189    !GG USE surface_data, ONLY : type_ocean, ok_veget
     190    USE surface_data, ONLY : type_ocean, ok_veget, iflag_seaice, amax_n, amax_s
     191    !GG
    167192    USE netcdf
    168193    USE indice_sol_mod
    … …  
    204229    REAL, DIMENSION(klon_mpi)                 :: rug_mpi  ! rugosity at global grid
    205230    REAL, DIMENSION(klon_mpi)                 :: alb_mpi  ! albedo at global grid
     231!GG
     232    REAL, DIMENSION(klon_glo)                 :: sih_glo  ! albedo at global grid
     233    REAL, DIMENSION(klon_mpi)                 :: sih_mpi  ! albedo at global grid
     234!GG
    206235
    207236    CHARACTER(len=20)                         :: modname='limit_read_mod'     
    … …  
    226255       END IF
    227256
     257       !GG
     258       IF (iflag_seaice==1) THEN
     259             ALLOCATE(sih(klon), stat=ierr)
     260             IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating sih',1)
     261       ENDIF
     262       !GG
     263
    228264       IF ( .NOT. ok_veget ) THEN
    229265          ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
    … …  
    304340         IF ( type_ocean /= 'couple') THEN                   
    305341             IF (is_omp_master) CALL xios_recv_field("sst_limin",sst_mpi)
     342             !GG
     343             IF (is_omp_master) CALL xios_recv_field("sih_limin",sih_mpi)
     344             !GG
    306345         ENDIF
    307346       
    … …  
    313352       IF ( type_ocean /= 'couple') THEN
    314353          CALL Scatter_omp(sst_mpi,sst)
     354          !GG
     355          CALL Scatter_omp(sih_mpi,sih)
     356          !GG
    315357          CALL Scatter_omp(pct_mpi(:,is_oce),pctsrf(:,is_oce))
    316358          CALL Scatter_omp(pct_mpi(:,is_sic),pctsrf(:,is_sic))
    … …  
    363405             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
    364406
     407! GG
     408! Account for leads
     409             IF (iflag_seaice>0) THEN
     410               DO ii=1,klon_glo/2
     411                 if (pct_glo(ii,is_sic)>amax_n) THEN
     412                    pct_glo(ii,is_oce)=pct_glo(ii,is_oce)+(pct_glo(ii,is_sic)-amax_n)
     413                    pct_glo(ii,is_sic)=amax_n
     414                 end if
     415               ENDDO
     416               DO ii=klon_glo/2,klon_glo
     417               if (pct_glo(ii,is_sic)>amax_s) THEN
     418                    pct_glo(ii,is_oce)=pct_glo(ii,is_oce)+(pct_glo(ii,is_sic)-amax_s)
     419                    pct_glo(ii,is_sic)=amax_s
     420               end if
     421               ENDDO
     422             ENDIF
     423!GG
    365424
    366425! Read land and continentals fraction only if asked for
    … …  
    395454             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
    396455             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
    397          
     456           !GG
     457             IF (iflag_seaice == 1) THEN
     458               ierr = NF90_INQ_VARID(nid, 'HICE', nvarid)
     459               IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <HICE> est absent',1)
     460
     461               ierr = NF90_GET_VAR(nid,nvarid,sih_glo(:),start,epais)
     462               IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <HICE>' ,1)
     463             ENDIF
     464            !GG
    398465          END IF
    399466
    … …  
    434501       IF ( type_ocean /= 'couple') THEN
    435502          CALL Scatter(sst_glo,sst)
     503          !GG
     504          IF (iflag_seaice==1) THEN
     505             CALL Scatter(sih_glo,sih)
     506          END IF
     507          !GG
    436508          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
    437509          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))

  • LMDZ6/trunk/libf/phylmd/ocean_forced_mod.F90

    r5486 r5662  
    2222       radsol, snow, agesno, &
    2323       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    24        tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa &
     24       tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa, &
     25       dthetadz300,pctsrf,Ampl &
    2526#ifdef ISO
    2627       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
    … …  
    3940    USE mod_grid_phy_lmdz
    4041    USE indice_sol_mod
     42    USE surface_data,     ONLY : iflag_leads
    4143    USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
    4244    use config_ocean_skin_m, only: activate_ocean_skin
    … …  
    6971    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
    7072    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
     73!GG
     74     REAL, DIMENSION(klon), INTENT(IN)        :: dthetadz300
     75     REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
     76!
    7177
    7278#ifdef ISO
    … …  
    9399    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
    94100    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
    95     REAL, intent(out):: sens_prec_liq(:) ! (knon)
     101    REAL, INTENT(out):: sens_prec_liq(:) ! (knon)
     102!GG
     103     REAL, DIMENSION(klon), INTENT(OUT)       :: Ampl
     104!
    96105
    97106#ifdef ISO     
    … …  
    110119    REAL, DIMENSION(knon)       :: sens_prec_sol
    111120    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol   
     121! GG
     122    REAL, DIMENSION(klon)       :: l_CBL, sicfra
     123!
    112124#ifdef ISO   
    113125    REAL, PARAMETER :: t_coup = 273.15     
    … …  
    204216    enddo
    205217
     218!GG
     219    if (iflag_leads == 1) then
     220      l_CBL = -52381.*dthetadz300 + 3008.1
     221      Ampl = 6.012e-08*l_CBL**2 - 4.036e-04*l_CBL + 1.4979
     222      WHERE(Ampl(:)>1.2) Ampl(:)=1.2
     223      sicfra(:)=pctsrf(:,is_sic)/(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
     224      WHERE(pctsrf(:,is_sic)+pctsrf(:,is_oce)<EPSFRA) sicfra(:)=0.
     225      WHERE(sicfra<0.7) Ampl(:)=1.
     226      WHERE((sicfra>0.7).and.(sicfra<0.9)) Ampl=((sicfra-0.7)/0.2)*Ampl+((0.9-sicfra)/0.2)
     227      fluxsens=Ampl*fluxsens
     228      dflux_s=Ampl*dflux_s
     229    endif
     230
    206231
    207232! - Flux calculation at first modele level for U and V
    … …  
    244269       AcoefH, AcoefQ, BcoefH, BcoefQ, &
    245270       AcoefU, AcoefV, BcoefU, BcoefV, &
    246        ps, u1, v1, gustiness, &
     271!GG       ps, u1, v1, gustiness, &
     272       ps, u1, v1, gustiness, pctsrf, &
     273!GG
    247274       radsol, snow, qsol, agesno, tsoil, &
    248275       qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    249        tsurf_new, dflux_s, dflux_l, rhoa &
     276!GG       tsurf_new, dflux_s, dflux_l, rhoa)
     277       tsurf_new, dflux_s, dflux_l, rhoa, swnet, hice, tice, bilg_cumul, &
     278       fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &
     279       dtice_melt, dtice_snow2sic &
     280!GG
    250281#ifdef ISO
    251282       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
    … …  
    261292    USE geometry_mod, ONLY: longitude,latitude
    262293    USE calcul_fluxs_mod
    263     USE surface_data,     ONLY : calice, calsno
     294!GG    USE surface_data,     ONLY : calice, calsno
     295    USE surface_data,     ONLY : calice, calsno, iflag_seaice, iflag_seaice_alb, &
     296            sice_cond, sisno_cond, sisno_den, sisno_min, sithick_min, sisno_wfact, &
     297            amax_s,amax_n, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, &
     298            si_pen_frac, si_pen_ext, fseaN, fseaS, iflag_leads
     299
     300    USE geometry_mod, ONLY: longitude,latitude,latitude_deg
     301!GG
    264302    USE limit_read_mod
    265303    USE fonte_neige_mod,  ONLY : fonte_neige
    … …  
    298336    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
    299337    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
     338!GG
     339    REAL, DIMENSION(klon), INTENT(IN)    :: swnet
     340    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
     341!GG
    300342#ifdef ISO
    301343    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
    … …  
    311353    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
    312354    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
     355!GG
     356    REAL, DIMENSION(klon), INTENT(INOUT)          :: hice
     357    REAL, DIMENSION(klon), INTENT(INOUT)          :: tice
     358    REAL, DIMENSION(klon), INTENT(INOUT)          :: bilg_cumul
     359    REAL, DIMENSION(klon), INTENT(INOUT)          :: fcds
     360    REAL, DIMENSION(klon), INTENT(INOUT)          :: fcdi
     361    REAL, DIMENSION(klon), INTENT(INOUT)          :: dh_basal_growth
     362    REAL, DIMENSION(klon), INTENT(INOUT)          :: dh_basal_melt
     363    REAL, DIMENSION(klon), INTENT(INOUT)          :: dh_top_melt
     364    REAL, DIMENSION(klon), INTENT(INOUT)          :: dh_snow2sic
     365    REAL, DIMENSION(klon), INTENT(INOUT)          :: dtice_melt
     366    REAL, DIMENSION(klon), INTENT(INOUT)          :: dtice_snow2sic
     367!GG
    313368#ifdef ISO     
    314369    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: xtsnow
    … …  
    342397    REAL, DIMENSION(knon)       :: sens_prec_liq, sens_prec_sol
    343398    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol   
     399!GG
     400    INTEGER                     :: ki
     401    INTEGER                     :: cpl_pas
     402    REAL, DIMENSION(klon)       :: bilg, fsic, f_bot
     403    REAL, PARAMETER             :: latent_ice = 334.0e3
     404    REAL, PARAMETER             :: rau_ice = 917.0
     405    REAL, PARAMETER             :: kice=2.2
     406    REAL                  :: f_cond, f_swpen, f_cond_s, f_cond_i
     407    REAL                  :: ustar, uscap, ustau
     408    ! for snow/ice albedo:
     409    REAL                  :: alb_snow, alb_ice, alb_pond
     410    REAL                  :: frac_snow, frac_ice, frac_pond
     411    REAL                  :: z1_i, z2_i, z1_s, zlog ! height parameters
     412    ! for ice melt / freeze
     413    REAL                  :: e_melt, snow_evap, h_test
     414    ! dhsic, dfsic change in ice mass, fraction.
     415    REAL                  :: dhsic, dfsic, frac_mf
     416    REAL                        :: fsea, amax
     417    REAL                  :: hice_i, tice_i, fsic_new
     418! snow and ice physical characteristics:
     419    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
     420    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
     421    REAL :: sno_den!=sisno_den !mean snow density, kg/m3
     422    REAL, PARAMETER :: ice_den=917. ! ice density
     423    REAL, PARAMETER :: sea_den=1025. ! sea water density
     424    REAL :: ice_cond!=sice_cond*ice_den !conductivity of ice
     425    REAL :: sno_cond!=sisno_cond*sno_den ! conductivity of snow
     426    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
     427    REAL, PARAMETER :: sea_cap=3995.   ! specific heat capacity, water
     428    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
     429
     430! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
     431    REAL :: snow_min!=sisno_min*sno_den !critical snow height 5 cm
     432    REAL :: snow_wfact!=sisno_wfact ! max fraction of falling snow blown into ocean
     433    REAL, PARAMETER :: ice_frac_min=0.005
     434    REAL :: h_ice_min!=sithick_min ! min ice thickness
     435    ! below ice_thin, priority is melt lateral / grow height
     436    ! ice_thin is also height of new ice
     437    REAL, PARAMETER :: h_ice_max=7 ! max ice height
     438    ! Ice thickness parameter for lateral growth
     439    REAL, PARAMETER :: h_ice_thick=1.5
     440    REAL, PARAMETER :: h_ice_thin=0.15
     441
     442! albedo  and radiation parameters
     443    INTEGER, SAVE :: iflag_sic_albedo
     444! albedo old or NEMO
     445    REAL :: alb_sno_dry!=rn_alb_sdry !dry snow albedo
     446    REAL :: alb_sno_wet!=rn_alb_smlt !wet snow albedo
     447    REAL :: alb_ice_dry!=rn_alb_idry !dry thick ice
     448    REAL :: alb_ice_wet!=rn_alb_imlt !melting thick ice
     449! new (Toyoda 2020) albedo
     450! Values for snow / ice, dry / melting, visible / near IR
     451    REAL, PARAMETER :: alb_sdry_vis=0.98
     452    REAL, PARAMETER :: alb_smlt_vis=0.88
     453    REAL, PARAMETER :: alb_sdry_nir=0.7
     454    REAL, PARAMETER :: alb_smlt_nir=0.55
     455    REAL, PARAMETER :: alb_idry_vis=0.78
     456    REAL, PARAMETER :: alb_imlt_vis=0.705
     457    REAL, PARAMETER :: alb_idry_nir=0.36
     458    REAL, PARAMETER :: alb_imlt_nir=0.285
     459    REAL, PARAMETER :: h_ice_alb=0.5*ice_den ! height for full ice albedo
     460    REAL, PARAMETER :: h_sno_alb=0.02*300 ! height for control of snow fraction
     461
     462    REAL :: pen_frac !=si_pen_frac !fraction of shortwave penetrating into the
     463    ! ice (not snow). Should be visible only, not NIR
     464    REAL :: pen_ext !=si_pen_ext !extinction length of penetrating shortwave (m-1)
     465
     466! HF from ocean below ice
     467!    REAL, PARAMETER :: fseaN=2.0 ! NH
     468!    REAL, PARAMETER :: fseaS=4.0 ! SH   
     469!GG
    344470
    345471#ifdef ISO
    … …  
    371497    tsurf_tmp(:) = tsurf_in(:)
    372498
     499!GG
     500    IF (iflag_seaice==0) THEN ! Old LMDZ sea ice surface
     501!GG
    373502! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal
    374503    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
    … …  
    387516       WHERE (snow > 0.0) cal = RCPD * calsno
    388517    ENDIF
     518
     519!GG
     520    ELSEIF (iflag_seaice==2) THEN
     521
     522    sno_den=sisno_den !mean snow density, kg/m3
     523    ice_cond=sice_cond*ice_den !conductivity of ice
     524    sno_cond=sisno_cond*sno_den ! conductivity of snow
     525    snow_min=sisno_min*sno_den !critical snow height 5 cm
     526    snow_wfact=sisno_wfact ! max fraction of falling snow blown into ocean
     527    h_ice_min=sithick_min ! min ice thickness
     528    alb_sno_dry=rn_alb_sdry !dry snow albedo
     529    alb_sno_wet=rn_alb_smlt !wet snow albedo
     530    alb_ice_dry=rn_alb_idry !dry thick ice
     531    alb_ice_wet=rn_alb_imlt !melting thick ice
     532    pen_frac=si_pen_frac !fraction of shortwave penetrating into the
     533    pen_ext=si_pen_ext !extinction length of penetrating shortwave (m-1)
     534
     535    bilg(:)=0.
     536    dif_grnd(:)=0.
     537    beta(:) = 1.
     538    fsic(:) = pctsrf(:,is_sic)
     539    cpl_pas =  NINT(86400./dtime * 1.0) ! une fois par jour
     540
     541    ! Surface, snow-ice and ice-ocean fluxes.
     542! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd)
     543
     544    !write(*,*) 'radsol 1',radsol(1:100)
     545    DO i=1,knon
     546        ki=knindex(i)
     547        IF (snow(i).GT.snow_min) THEN
     548            !  1 / snow-layer heat capacity
     549            cal(i)=2.*RCPD/(snow(i)*ice_cap)
     550            ! adjustment time-scale of conductive flux
     551            dif_grnd(i) = cal(i) * sno_cond / snow(i) / RCPD
     552            ! for conductive flux
     553            f_cond_s = sno_cond * (tice(ki)-t_freeze) / snow(i)
     554            radsol(i) = radsol(i)+f_cond_s
     555            ! all shortwave flux absorbed
     556            f_swpen=0.
     557        ELSE ! bare ice.
     558            f_cond_s = 0.
     559            ! 1 / ice-layer heat capacity
     560            cal(i) = 2.*RCPD/(hice(ki)*ice_den*ice_cap)
     561            ! adjustment time-scale of conductive flux
     562            dif_grnd(i) = cal(i) * ice_cond / (hice(ki)*ice_den) / RCPD
     563            ! penetrative shortwave flux...
     564            f_swpen=swnet(i)*pen_frac*exp(-pen_ext*hice(ki))
     565            radsol(i) = radsol(i)-f_swpen
     566            ! GG no conductive flux in this case?
     567        END IF
     568        bilg(ki)=f_swpen-f_cond_s
     569    END DO
     570
     571    endif
    389572
    390573!    beta = 1.0
    … …  
    423606!
    424607!****************************************************************************************
     608!GG
     609    if (iflag_seaice==0) then
     610!GG
    425611#ifdef ISO
    426612   ! verif
    … …  
    502688    alb2_new(:) = alb1_new(:)
    503689
     690!GG
     691  else
     692
     693        DO i=1,knon
     694        ki=knindex(i)
     695
     696           ! ocean-ice heat flux
     697           fsea=fseaS
     698           amax=amax_s
     699           if (latitude(ki)>0) THEN
     700                   fsea=fseaN
     701                   amax=amax_n
     702           ENDIF
     703
     704           IF (snow(i).GT.snow_min) THEN
     705                ! snow conductive flux after calcul_fluxs (pos up)
     706                f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i)
     707                ! 1 / heat capacity and conductive timescale
     708                uscap = 2. / ice_cap / (snow(i)+hice(ki)*ice_den)
     709                ustau = uscap * ice_cond / (hice(ki)*ice_den)
     710                ! update ice temp
     711                tice(ki) = (tice(ki) + dtime*(ustau*t_freeze - uscap*f_cond_s)) &
     712                     / (1 + dtime*ustau)
     713           ELSE ! bare ice
     714                tice(ki)=tsurf_new(i)
     715           ENDIF
     716           ! ice conductive flux (pos up)
     717           f_cond_i = ice_cond * (t_freeze-tice(ki)) / (hice(ki)*ice_den)
     718           f_bot(i) = fsea - f_cond_i
     719           fcdi(ki) = f_cond_i - fsea
     720           fcds(i) = f_cond_s
     721           !bilg(ki) = bilg(ki)+f_cond_i
     722        END DO
     723
     724!****************************************************************************************
     725! 2) Update snow and ice surface : thickness
     726!****************************************************************************************
     727
     728    IF (iflag_seaice==1) THEN
     729!   Read from limit
     730    CALL limit_read_hice(knon,knindex,hice)
     731    ENDIF
     732!   Formula Krinner et al. 1997 : h = (0.2 + 3.8(f_min**2))(1 + 2(f- f_min))
     733
     734
     735
     736    DO i=1,knon
     737        ki=knindex(i)
     738        IF (precip_snow(i) > 0.) THEN
     739            snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
     740        END IF
     741! snow and ice sublimation
     742        IF (evap(i) > 0.) THEN
     743           snow_evap = MIN (snow(i) / dtime, evap(i))
     744           snow(i) = snow(i) - snow_evap * dtime
     745           snow(i) = MAX(0.0, snow(i))
     746           IF (iflag_seaice==2) THEN
     747             hice(ki) = MAX(0.0,hice(ki)-(evap(i)-snow_evap)*dtime/ice_den)
     748           ENDIF
     749        ENDIF
     750! Melt / Freeze snow from above if Tsurf>0
     751        IF (tsurf_new(i).GT.t_melt) THEN
     752            ! energy available for melting snow (in kg of melted snow /m2)
     753            e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
     754               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
     755            ! remove snow
     756            tice_i=tice(ki)
     757            IF (snow(i).GT.e_melt) THEN
     758                snow(i)=snow(i)-e_melt
     759                tsurf_new(i)=t_melt
     760            ELSE ! all snow is melted
     761                ! add remaining heat flux to ice
     762                e_melt=e_melt-snow(i)
     763                tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*hice(ki)*ice_den)
     764                tsurf_new(i)=tice(ki)
     765            END IF
     766            dtice_melt(ki)=(tice(ki)-tice_i)/dtime
     767        END IF
     768! Bottom melt / grow
     769! bottom freeze if bottom flux (cond + oce-ice) <0
     770        IF (iflag_seaice==2) THEN
     771         IF (f_bot(i).LT.0) THEN
     772           ! larger fraction of bottom growth
     773           frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thick)   &
     774                  / (h_ice_max-h_ice_thick)))
     775           ! quantity of new ice (formed at mean ice temp)
     776           e_melt= -f_bot(i) * dtime * fsic(ki) &
     777                   / (ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     778           ! first increase height to h_thick
     779           dhsic=MAX(0.,MIN(h_ice_thick-hice(ki),e_melt/(fsic(ki)*ice_den)))
     780           hice_i=hice(ki)
     781           hice(ki)=dhsic+hice(ki)
     782           e_melt=e_melt-fsic(ki)*dhsic
     783           IF (e_melt.GT.0.) THEN
     784           ! frac_mf fraction used for lateral increase
     785           dfsic=MIN(amax-fsic(ki),e_melt*frac_mf/ (hice(ki)*ice_den) )
     786           ! No lateral growth -> forced ocean
     787           !fsic(ki)=fsic(ki)+dfsic
     788           e_melt=e_melt-dfsic*hice(ki)*ice_den
     789           ! rest used to increase height
     790           hice(ki)=MIN(h_ice_max,hice(ki)+e_melt/( fsic(ki) * ice_den ) )
     791           END IF
     792           dh_basal_growth(ki)=(hice(ki)-hice_i)/dtime
     793
     794! melt from below if bottom flux >0
     795         ELSE
     796           ! larger fraction of lateral melt from warm ocean
     797           frac_mf=MIN(1.,MAX(0.,(hice(ki)-h_ice_thin)   &
     798                  / (h_ice_thick-h_ice_thin)))
     799           ! bring ice to freezing and melt from below
     800           ! quantity of melted ice
     801           e_melt= f_bot(i) * dtime * fsic(ki) &
     802                   / (ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
     803           ! first decrease height to h_thick
     804           hice_i=hice(ki)
     805           dhsic=MAX(0.,MIN(hice(ki)-h_ice_thick,e_melt/(fsic(ki)*ice_den)))
     806           hice(ki)=hice(ki)-dhsic
     807           e_melt=e_melt-fsic(ki)*dhsic*ice_den
     808
     809           IF (e_melt.GT.0) THEN
     810           ! frac_mf fraction used for height decrease
     811           dhsic=MAX(0.,MIN(hice(ki)-h_ice_min,e_melt/ice_den*frac_mf/fsic(ki)))
     812           hice(ki)=hice(ki)-dhsic
     813           e_melt=e_melt-fsic(ki)*dhsic*ice_den
     814           ! rest used to decrease fraction (up to 0!)
     815           dfsic=MIN(fsic(ki),e_melt/(hice(ki)*ice_den))
     816           ! Remaining heat not used if everything melted
     817           e_melt=e_melt-dfsic*hice(ki)*ice_den
     818           ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
     819           END IF
     820           dh_basal_melt(ki)=(hice(ki)-hice_i)/dtime
     821         END IF
     822        END IF
     823
     824! melt ice from above if Tice>0
     825        tice_i=tice(ki)
     826        IF (tice(ki).GT.t_melt) THEN
     827           IF (iflag_seaice==2) THEN
     828            ! quantity of ice melted (kg/m2)
     829            e_melt=MAX(hice(ki)*ice_den*(tice(ki)-t_melt)*ice_cap/2. &
     830             /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
     831            ! melt from above, height only
     832            hice_i=hice(ki)
     833            dhsic=MIN(hice(ki)-h_ice_min,e_melt/ice_den)
     834            dh_top_melt(i)=dhsic
     835            e_melt=e_melt-dhsic
     836            IF (e_melt.GT.0) THEN
     837              ! lateral melt if ice too thin
     838              dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/(h_ice_min*ice_den)*fsic(ki))
     839              ! if all melted do nothing with remaining heat
     840              e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min*ice_den)
     841              ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
     842            END IF
     843            hice(ki)=hice(ki)-dhsic
     844            dh_top_melt(ki)=(hice(ki)-hice_i)/dtime
     845            ! surface temperature at melting point
     846           END IF
     847           tice(ki)=t_melt
     848           tsurf_new(i)=t_melt
     849        END IF
     850        dtice_melt(ki)=dtice_melt(ki)+tice(ki)-tice_i
     851
     852        ! convert snow to ice if below floating line
     853        h_test=(hice(ki)*ice_den+snow(i))-hice(ki)*sea_den
     854        IF ((h_test.GT.0.).AND.(hice(ki).GT.h_ice_min)) THEN !snow under water
     855            ! extra snow converted to ice (with added frozen sea water)
     856            IF (iflag_seaice==2) THEN
     857             dhsic=h_test/(sea_den-ice_den+sno_den)
     858             hice(ki)=hice(ki)+dhsic
     859            ENDIF
     860            snow(i)=snow(i)-dhsic*sno_den
     861            ! available energy (freeze sea water + bring to tice)
     862            e_melt=dhsic*ice_den*(1.-sno_den/ice_den)*(ice_lat+ &
     863                   ice_cap/2.*(t_freeze-tice(ki)))
     864            ! update ice temperature
     865            tice_i=tice(ki)
     866            tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+hice(ki)*ice_den)
     867            IF (iflag_seaice==2) THEN
     868              dh_snow2sic(ki)=dhsic/dtime
     869            END IF
     870            dtice_snow2sic(ki)=(tice(ki)-tice_i)/dtime
     871        END IF
     872    END DO
     873
     874        !write(*,*) 'hice 2',hice(1:100)
     875        !write(*,*) 'tice 2',tice(1:100)
     876
     877        iflag_sic_albedo=iflag_seaice_alb
     878
     879!*******************************************************************************
     880! 3) cumulate ice-ocean fluxes, update tslab, lateral grow
     881!***********************************************o*******************************
     882    !cumul fluxes.
     883    bilg_cumul(:)=bilg_cumul(:)+bilg(:)/float(cpl_pas)
     884    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab
     885        bilg_cumul(:)=0.
     886    END IF ! coupling time
     887
     888!   write(*,*) 'hice 3',hice(1:100)
     889!   write(*,*) 'tice 3',tice(1:100)
     890    !tests ice fraction
     891    WHERE (fsic.LT.ice_frac_min)
     892        tice=t_melt
     893        hice=h_ice_min
     894    END WHERE
     895
     896    !write(*,*) 'hice 4',hice(1:100)
     897    !write(*,*) 'tice 4',tice(1:100)
     898
     899    endif
     900
     901!****************************************************************************************
     902! 4) Compute sea-ice and snow albedo
     903!****************************************************************************************
     904    IF (iflag_seaice > 0) THEN
     905    SELECT CASE (iflag_sic_albedo)
     906      CASE(0)
     907! old slab albedo : single value. age of snow, melt ponds.
     908        DO i=1,knon
     909          ki=knindex(i)
     910         ! snow albedo: update snow age
     911          IF (snow(i).GT.0.0001) THEN
     912               agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
     913                           * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.)
     914          ELSE
     915              agesno(i)=0.0
     916          END IF
     917          ! snow albedo
     918          alb_snow=alb_sno_wet+(alb_sno_dry-alb_sno_wet)*EXP(-agesno(i)/50.)
     919          ! ice albedo (varies with ice tkickness and temp)
     920          alb_ice=MAX(0.0,0.13*LOG(100.*hice(ki))+0.1)
     921          !alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1)
     922          IF (tice(ki).GT.t_freeze-0.01) THEN
     923              alb_ice=MIN(alb_ice,alb_ice_wet)
     924          ELSE
     925              alb_ice=MIN(alb_ice,alb_ice_dry)
     926          END IF
     927          ! pond albedo
     928          alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
     929          ! pond fraction
     930          frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
     931          ! snow fraction
     932          frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
     933          ! ice fraction
     934          frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
     935          ! total albedo
     936          alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond
     937        END DO
     938        alb2_new(:) = alb1_new(:)
     939
     940      CASE(1)
     941! New visible and IR albedos, dry / melting snow
     942! based on Toyoda et al, 2020
     943      DO i=1,knon
     944          ki=knindex(i)
     945          ! snow fraction
     946          frac_snow  = snow(i) / (snow(i) + h_sno_alb)
     947          ! dependence of ice albedo with ice thickness
     948          frac_ice = MIN(1.,ATAN(4.*hice(ki)*ice_den) / ATAN(4.*h_ice_alb))
     949          ! Total (for ice, min = 0.066 = alb_ocean)
     950          IF (tice(ki).GT.t_melt) THEN
     951              alb_ice = 0.066 + (alb_imlt_vis - 0.066)*frac_ice
     952              alb1_new(i)=alb_smlt_vis*frac_snow + alb_ice*(1.-frac_snow)
     953              alb_ice = 0.066 + (alb_imlt_nir - 0.066)*frac_ice
     954              alb2_new(i)=alb_smlt_nir*frac_snow + alb_ice*(1.-frac_snow)
     955          ELSEIF (tice(ki).GT.t_melt - 1.) THEN
     956              frac_pond = tice(ki) - t_freeze
     957              alb_snow = alb_smlt_vis*frac_pond + alb_sdry_vis*(1.-frac_pond)
     958              alb_ice = alb_imlt_vis*frac_pond + alb_idry_vis*(1.-frac_pond)
     959              alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice
     960              alb1_new(i)= alb_snow*frac_snow + alb_ice*(1.-frac_snow)
     961              alb_snow = alb_smlt_nir*frac_pond + alb_sdry_nir*(1.-frac_pond)
     962              alb_ice = alb_imlt_nir*frac_pond + alb_idry_nir*(1.-frac_pond)
     963              alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice
     964              alb2_new(i)= alb_snow*frac_snow + alb_ice*(1.-frac_snow)
     965          ELSE
     966              alb_ice = 0.066 + (alb_idry_vis - 0.066)*frac_ice
     967              alb1_new(i)=alb_sdry_vis*frac_snow + alb_ice*(1.-frac_snow)
     968              alb_ice = 0.066 + (alb_idry_nir - 0.066)*frac_ice
     969              alb2_new(i)=alb_sdry_nir*frac_snow + alb_ice*(1.-frac_snow)
     970          ENDIF
     971      END DO
     972
     973      CASE(2)
     974! LIM3 scheme. Uses clear sky / overcast value, with 50% clear sky
     975      z1_i = 1.5 * ice_den
     976      z2_i = 0.05 * ice_den
     977      zlog = 1. / (LOG(1.5) - LOG(0.05))
     978      z1_s = 1. / (0.025 * sno_den)
     979      DO i=1,knon
     980          ki=knindex(i)
     981            ! temperature above / below 0
     982            IF (tice(ki).GE.t_melt) THEN
     983               alb_ice = alb_ice_wet
     984               alb_snow = alb_sno_wet
     985            ELSE
     986               alb_ice = alb_ice_dry
     987               alb_snow = alb_sno_dry
     988            ENDIF
     989            ! ice thickness
     990            IF (hice(ki)*ice_den.LT.z2_i) THEN
     991                alb_ice = 0.066 + 0.114 * hice(ki)*ice_den / z2_i
     992            ELSEIF (hice(ki)*ice_den.LT.z1_i) THEN
     993                alb_ice = alb_ice + (0.18 - alb_ice) &
     994                          * (LOG(z1_i) - LOG(hice(ki)*ice_den)) * zlog
     995            ENDIF
     996            ! ice or snow depending on snow thickness
     997            alb_snow = alb_snow - (alb_snow -alb_ice) * EXP(- snow(i) * z1_s)
     998            ! Effect of clouds (polynomial fit with 50% clouds)
     999            alb1_new(i) = alb_snow - 0.5 * (-0.1010 * alb_snow*alb_snow &
     1000                          + 0.1933*alb_snow - 0.0148)
     1001            alb2_new(i) = alb1_new(i)
     1002      END DO
     1003
     1004      CASE(3)
     1005      CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))
     1006      WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.
     1007      alb1_new(:) = 0.0
     1008      DO i=1, knon
     1009         zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
     1010         alb1_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
     1011      ENDDO
     1012      alb2_new(:) = alb1_new(:)
     1013
     1014      print*,'alb_neig=',alb_neig
     1015      print*,'zfra=',zfra
     1016      print*,'snow=',snow
     1017      print*,'alb1_new=',alb1_new
     1018      print*,'alb2_new=',alb2_new
     1019    END SELECT
     1020    END IF
     1021! ------ End Albedo ----------
     1022
     1023!GG
    5041024  END SUBROUTINE ocean_forced_ice
    5051025

  • LMDZ6/trunk/libf/phylmd/ocean_slab_mod.f90

    r5285 r5662  
    5050  !$OMP THREADPRIVATE(fsic)
    5151  ! temperature of the sea ice
    52   REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
    53   !$OMP THREADPRIVATE(tice)
     52!GG
     53  !REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice_slab
     54  !!$OMP THREADPRIVATE(tice_slab)
     55  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice_slab
     56  !$OMP THREADPRIVATE(tice_slab)
     57!GG
    5458  ! sea ice thickness, in kg/m2
    5559  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
    … …  
    238242        ENDIF
    239243        bilg_cum(:) = 0.0   
    240         ALLOCATE(tice(klon), stat = error)
     244!GG        ALLOCATE(tice(klon), stat = error)
     245        ALLOCATE(tice_slab(klon), stat = error)
    241246        IF (error /= 0) THEN
    242247           abort_message='Pb allocation slab_tice'
    … …  
    633638                      e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
    634639                      ! new ice
    635                       tice(ki)=t_freeze
     640                !GG      tice(ki)=t_freeze
     641                      tice_slab(ki)=t_freeze
    636642                      fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
    637643                      IF (fsic(ki).GT.ice_frac_min) THEN
    … …  
    650656                      ! quantity of new ice formed over open ocean
    651657                      e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
    652                                /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     658                               /(ice_lat+ice_cap/2.*(t_freeze-tice_slab(ki)))
     659              !GG                 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
    653660                      ! new ice height and fraction
    654661                      h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
    … …  
    759766            cal(i)=2.*RCPD/(snow(i)*ice_cap)
    760767            ! snow conductive flux
    761             f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
     768            f_cond=sno_cond*(tice_slab(ki)-tsurf_in(i))/snow(i)
     769       !GG     f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
    762770            ! all shortwave flux absorbed
    763771            f_swpen=0.
    764772            ! bottom flux (ice conduction)
    765             slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
     773            slab_bilg(ki)=ice_cond*(tice_slab(ki)-t_freeze)/seaice(ki)
     774       !GG     slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
    766775            ! update ice temperature
    767             tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
     776       !GG     tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
     777            tice_slab(ki)=tice_slab(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
    768778                     *(slab_bilg(ki)+f_cond)*dtime
    769779       ELSE ! bare ice
    … …  
    771781            cal(i)=2.*RCPD/(seaice(ki)*ice_cap)
    772782            ! conductive flux
    773             f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
     783       !GG     f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
     784            f_cond=ice_cond*(t_freeze-tice_slab(ki))/seaice(ki)
    774785            ! penetrative shortwave flux...
    775786            f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den)
    … …  
    789800        tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    790801    DO i=1,knon
    791         IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
     802    !GG    IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
     803        IF (snow(i).LT.snow_min) tice_slab(knindex(i))=tsurf_new(i)
    792804    END DO
    793805
    … …  
    819831            ! energy available for melting snow (in kg of melted snow /m2)
    820832            e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
    821                /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
     833               /(ice_lat+ice_cap/2.*(t_melt-tice_slab(ki))),0.0),snow(i))
     834        !GG       /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
    822835            ! remove snow
    823836            IF (snow(i).GT.e_melt) THEN
    … …  
    827840                ! add remaining heat flux to ice
    828841                e_melt=e_melt-snow(i)
    829                 tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
    830                 tsurf_new(i)=tice(ki)
     842                tice_slab(ki)=tice_slab(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
     843         !GG       tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
     844                tsurf_new(i)=tice_slab(ki)
     845         !GG       tsurf_new(i)=tice(ki)
    831846            END IF
    832847        END IF
    833848! melt ice from above if Tice>0
    834         IF (tice(ki).GT.t_melt) THEN
     849      !GG  IF (tice(ki).GT.t_melt) THEN
     850        IF (tice_slab(ki).GT.t_melt) THEN
    835851            ! quantity of ice melted (kg/m2)
    836             e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. &
     852      !GG      e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. &
     853            e_melt=MAX(seaice(ki)*(tice_slab(ki)-t_melt)*ice_cap/2. &
    837854             /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
    838855            ! melt from above, height only
    … …  
    850867            seaice(ki)=seaice(ki)-dhsic
    851868            ! surface temperature at melting point
    852             tice(ki)=t_melt
     869        !GG    tice(ki)=t_melt
     870            tice_slab(ki)=t_melt
    853871            tsurf_new(i)=t_melt
    854872        END IF
    … …  
    860878            seaice(ki)=seaice(ki)+dhsic
    861879            snow(i)=snow(i)-dhsic*sno_den/ice_den
    862             ! available energy (freeze sea water + bring to tice)
     880            ! available energy (freeze sea water + bring to tice_slab)
    863881            e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
    864                    ice_cap/2.*(t_freeze-tice(ki)))
     882                   ice_cap/2.*(t_freeze-tice_slab(ki)))
     883       !GG            ice_cap/2.*(t_freeze-tice(ki)))
    865884            ! update ice temperature
    866             tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
     885       !GG     tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
     886            tice_slab(ki)=tice_slab(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
    867887        END IF
    868888    END DO
    … …  
    882902        ! ice albedo (varies with ice tkickness and temp)
    883903        alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1)
    884         IF (tice(ki).GT.t_freeze-0.01) THEN
     904     !GG   IF (tice(ki).GT.t_freeze-0.01) THEN
     905        IF (tice_slab(ki).GT.t_freeze-0.01) THEN
    885906            alb_ice=MIN(alb_ice,alb_ice_wet)
    886907        ELSE
    … …  
    888909        END IF
    889910        ! pond albedo
    890         alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
     911        alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice_slab(ki)-t_melt,-2.0)))
     912     !GG   alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
    891913        ! pond fraction
    892         frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
     914        frac_pond=0.2*(2.0+MIN(0.0,MAX(tice_slab(ki)-t_melt,-2.0)))
     915     !GG   frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
    893916        ! snow fraction
    894917        frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
    … …  
    917940               ! quantity of new ice
    918941                e_melt=(t_freeze-tslab(ki,1))/cyang &
    919                        /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     942                       /(ice_lat+ice_cap/2.*(t_freeze-tice_slab(ki)))
     943       !GG                /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
    920944               ! first increase height to h_thin
    921945               dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
    … …  
    935959               ! quantity of melted ice
    936960               e_melt=(tslab(ki,1)-t_freeze)/cyang &
    937                        /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
     961                       /(ice_lat+ice_cap/2.*(tice_slab(ki)-t_freeze))
     962          !GG             /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
    938963               ! first decrease height to h_thick
    939964               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
    … …  
    960985    WHERE (fsic.LT.ice_frac_min)
    961986        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
    962         tice=t_melt
     987        tice_slab=t_melt
     988   !GG     tice=t_melt
    963989        fsic=0.
    964990        seaice=0.
    … …  
    9761002    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
    9771003    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
    978     IF (ALLOCATED(tice)) DEALLOCATE(tice)
     1004    IF (ALLOCATED(tice_slab)) DEALLOCATE(tice_slab)
    9791005    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
    9801006    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)

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

    r5627 r5662  
    1414  USE mod_grid_phy_lmdz,   ONLY : klon_glo
    1515  USE ioipsl
    16   USE surface_data,        ONLY : type_ocean, ok_veget, landice_opt
     16  USE surface_data,        ONLY : type_ocean, ok_veget, landice_opt, iflag_leads
    1717  USE surf_land_mod,       ONLY : surf_land
    1818  USE surf_landice_mod,    ONLY : surf_landice
    … …  
    4242  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
    4343  !$OMP THREADPRIVATE(fder)
     44!GG
     45  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: hice   ! flux drift
     46  !$OMP THREADPRIVATE(hice)
     47  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: tice   ! flux drift
     48  !$OMP THREADPRIVATE(tice)
     49  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: bilg_cumul   ! flux drift
     50  !$OMP THREADPRIVATE(bilg_cumul)
     51!GG
    4452  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE    :: snow   ! snow at surface
    4553  !$OMP THREADPRIVATE(snow)
    … …  
    7684!****************************************************************************************
    7785!
    78   SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
     86!GG
     87!  SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
     88  SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst, hice_rst,tice_rst,bilg_cumul_rst)
     89!GG
    7990
    8091! This routine should be called after the restart file has been read.
    … …  
    91102!****************************************************************************************
    92103    REAL, DIMENSION(klon), INTENT(IN)                 :: fder_rst
     104!GG
     105    REAL, DIMENSION(klon), INTENT(IN)                 :: hice_rst
     106    REAL, DIMENSION(klon), INTENT(IN)                 :: tice_rst
     107    REAL, DIMENSION(klon), INTENT(IN)                 :: bilg_cumul_rst
     108!GG
    93109    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: snow_rst
    94110    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: qsurf_rst
    … …  
    108124    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
    109125
     126!GG
     127    ALLOCATE(hice(klon), stat=ierr)
     128    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
     129
     130    ALLOCATE(tice(klon), stat=ierr)
     131    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
     132
     133    ALLOCATE(bilg_cumul(klon), stat=ierr)
     134    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
     135!GG
     136
    110137    ALLOCATE(snow(klon,nbsrf), stat=ierr)
    111138    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
    … …  
    124151
    125152    fder(:)       = fder_rst(:)
     153!GG
     154    hice(:)       = hice_rst(:)
     155    tice(:)       = tice_rst(:)
     156    bilg_cumul(:)       = bilg_cumul_rst(:)
     157!GG
    126158    snow(:,:)     = snow_rst(:,:)
    127159    qsurf(:,:)    = qsurf_rst(:,:)
    … …  
    261293       debut,     lafin,                              &
    262294       rlon,      rlat,      rugoro,   rmu0,          &
    263        lwdown_m,  cldt,                               &
     295   !GG lwdown_m,  cldt,          &
     296       lwdown_m,  pphi, cldt,          &
     297   !GG
    264298       rain_f,    snow_f,    bs_f, solsw_m,  solswfdiff_m, sollw_m,       &
    265299       gustiness,                                     &
    … …  
    312346!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
    313347!!        tke_x,     tke_w                              &
    314        wake_dltke,                                   &
    315         treedrg,                                      &
     348       wake_dltke,                                     &
     349!GG        treedrg                                   &
     350       treedrg,hice ,tice, bilg_cumul,            &
     351       fcds, fcdi, dh_basal_growth, dh_basal_melt, &
     352       dh_top_melt, dh_snow2sic, &
     353       dtice_melt, dtice_snow2sic , &
     354!GG
    316355!FC
    317356!AM heterogeneous continental sub-surfaces
    … …  
    369408! cldt-----input-R- total cloud fraction
    370409! Martin
     410!GG
     411! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol)
     412!GG
    371413!
    372414! d_t------output-R- le changement pour "t"
    … …  
    470512    REAL, DIMENSION(klon),        INTENT(IN)        :: gustiness ! gustiness
    471513
     514!GG
     515    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pphi    ! geopotential (m2/s2)
     516!GG
    472517    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud
    473518
    … …  
    679724    REAL, DIMENSION(klon),       INTENT(OUT)        :: runoff     ! runoff on land ice
    680725! Martin
     726!GG
     727    REAL, DIMENSION(klon),       INTENT(INOUT)        :: hice      ! hice
     728    REAL, DIMENSION(klon),       INTENT(INOUT)        :: tice      ! tice
     729    REAL, DIMENSION(klon),       INTENT(INOUT)        :: bilg_cumul      ! flux cumulated
     730    REAL, DIMENSION(klon),       INTENT(INOUT)        :: fcds
     731    REAL, DIMENSION(klon),       INTENT(INOUT)        :: fcdi
     732    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dh_basal_growth
     733    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dh_basal_melt
     734    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dh_top_melt
     735    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dh_snow2sic
     736    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dtice_melt
     737    REAL, DIMENSION(klon),       INTENT(INOUT)        :: dtice_snow2sic
     738!GG
    681739
    682740! Local variables with attribute SAVE
    … …  
    10751133    ! dt_ds, tkt, tks, taur, sss on ocean points
    10761134    REAL :: missing_val
     1135
     1136    ! GG
     1137    REAL, DIMENSION(klon,klev)         :: ytheta
     1138    REAL, DIMENSION(klon,klev)         :: ypphii
     1139    REAL, DIMENSION(klon,klev)         :: ypphi
     1140    REAL, DIMENSION(klon,klev)         :: ydthetadz
     1141    REAL, DIMENSION(klon)              :: ydthetadz300
     1142    REAL, DIMENSION(klon)              :: Ampl
     1143    ! GG
     1144
    10771145    ! AM !
    10781146    REAL, DIMENSION(klon) :: z0m_eff, z0h_eff, ratio_z0m_z0h_eff, albedo_eff
    … …  
    13981466   yfields_out(:,:) = 0.
    13991467! << PC
     1468
     1469!GG
     1470  ypphi = 0.0 
     1471!GG
    14001472
    14011473
    … …  
    17961868             yq(j,k) = q(i,k)
    17971869             yqbs(j,k)=qbs(i,k)
     1870!! GG
     1871             ypphi(j,k) = pphi(i,k)
     1872!!
     1873
    17981874#ifdef ISO
    17991875             DO ixt=1,ntraciso   
    … …  
    24912567               cdragm_tersrf, cdragh_tersrf, &
    24922568               swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf  &
     2569!GG
     2570!               yveget,ylai,yheight,hice,tice,bilg_cumul, &
     2571!               fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &
     2572!               dtice_melt, dtice_snow2sic)
     2573               !GG
    24932574#ifdef ISO
    24942575         &      ,yxtrain_f, yxtsnow_f,yxt1, &
    … …  
    26252706         
    26262707       CASE(is_oce)
     2708
     2709!GG
     2710! calculate length scale PBL
     2711
     2712        if (iflag_leads == 1) then
     2713        ydthetadz = 999999.
     2714        ypphii = 999999.
     2715        ytheta = 999999.
     2716
     2717        DO k = 1, klev
     2718          DO j = 1, knon
     2719             ytheta(j,k) = yt(j,k)*(ypplay(j,k)/1.e5)**(RD/RCPD)
     2720          ENDDO
     2721        ENDDO
     2722
     2723        DO k = 2, klev
     2724          DO j = 1, knon
     2725             ydthetadz(j,k) = RG*( ytheta(j,k) - ytheta(j,k-1) ) / ( ypphi(j,k) - ypphi(j,k-1) )
     2726             ypphii(j,k) = (ypphi(j,k)+ypphi(j,k-1))/(RG*2.)
     2727          ENDDO
     2728        ENDDO
     2729
     2730        DO j = 1, knon
     2731            ! print *, "ypphii(j,:)=", ypphii(j,:)
     2732            ! print *, "ypplay(j,:)=", ypplay(j,:)
     2733            ! print *, "ytheta(j,:)=", ytheta(j,:)
     2734            ! print *, "minloc(abs(ypphii(j,:)-300))=",
     2735            ! minloc(abs(ypphii(j,:)-300),1)
     2736             k= minloc(abs(ypphii(j,:)-300),1)
     2737             ydthetadz300(j)=ydthetadz(j,k)
     2738        ENDDO
     2739        end if
     2740!GG
    26272741           CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
    26282742               ywindsp, rmu0, yfder, yts, &
    … …  
    26382752               y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
    26392753               yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
    2640                ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss &
     2754           !GG    ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss)
     2755               ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss, &
     2756               ydthetadz300,Ampl                 &
     2757           !GG
    26412758#ifdef ISO
    26422759         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
    … …  
    26842801!albedo SB <<<
    26852802               ytsurf_new, y_dflux_t, y_dflux_q, &
    2686                y_flux_u1, y_flux_v1 &
     2803!GG               y_flux_u1, y_flux_v1)
     2804               y_flux_u1, y_flux_v1, &
     2805               hice,tice,bilg_cumul, &
     2806               fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &
     2807               dtice_melt, dtice_snow2sic     &
     2808!GG
    26872809#ifdef ISO
    26882810         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &

  • LMDZ6/trunk/libf/phylmd/phyaqua_mod.f90

    r5645 r5662  
    330330
    331331
    332     CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
     332!GG
     333    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil, hice, tice, bilg_cumul)
     334!    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
     335!GG
    333336
    334337    PRINT *, 'iniaqua: before phyredem'

  • LMDZ6/trunk/libf/phylmd/phyetat0_mod.f90

    r5645 r5662  
    1616  USE fonte_neige_mod,  ONLY : fonte_neige_init
    1717  USE pbl_surface_mod,  ONLY : pbl_surface_init
    18   USE surface_data,     ONLY : type_ocean, version_ocean
     18!GG  USE surface_data,     ONLY : type_ocean, version_ocean
     19  USE surface_data,     ONLY : type_ocean, version_ocean, iflag_seaice, &
     20                                   iflag_seaice_alb, iflag_leads
     21!GG
    1922  USE phyetat0_get_mod, ONLY : phyetat0_get, phyetat0_srf
    2023  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, phys_tstep, &
    … …  
    3033       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
    3134       ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, dter, dser, &
    32        dt_ds, ratqs_inter_, frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, &
     35!GG       dt_ds, ratqs_inter_
     36       dt_ds, ratqs_inter_, &
     37       hice, tice, bilg_cumul, &
     38!GG
     39       frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, &
    3340       albedo_tersrf, beta_tersrf, inertie_tersrf, alpha_soil_tersrf, &
    3441       period_tersrf, hcond_tersrf, tsurfi_tersrf, tsoili_tersrf, tsoil_depth, &
    … …  
    4350  USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl, carbon_cycle_tr, carbon_cycle_rad, co2_send, RCO2_glo
    4451  USE indice_sol_mod,   ONLY: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
    45   USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice, ocean_slab_init
     52  !GG USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice, ocean_slab_init
     53  USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice_slab, ocean_slab_init
     54  !GG
    4655  USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
    4756  use wxios_mod, ONLY: missing_val_xios => missing_val, using_xios
    … …  
    162171  IF (ok_orolf) tab_cntrl(11) =1.
    163172  IF (ok_limitvrai) tab_cntrl(12) =1.
     173  !GG
     174  tab_cntrl(18) =iflag_seaice
     175  tab_cntrl(19) =iflag_seaice_alb
     176  tab_cntrl(20) =iflag_leads
     177  !GG
    164178
    165179  itau_phy = tab_cntrl(15)
    … …  
    637651      ! Sea ice variables
    638652      IF (version_ocean == 'sicINT') THEN
    639           found=phyetat0_get(tice,"slab_tice","slab_tice",0.)
     653          found=phyetat0_get(tice_slab,"slab_tice","slab_tice",0.)
     654  !GG        found=phyetat0_get(tice,"slab_tice","slab_tice",0.)
    640655          IF (.NOT. found) THEN
    641               PRINT*, "phyetat0: Le champ <tice> est absent"
     656  !GG            PRINT*, "phyetat0: Le champ <tice> est absent"
     657              PRINT*, "phyetat0: Le champ <tice_slab> est absent"
    642658              PRINT*, "Initialisation a tsol_sic"
    643                   tice(:)=ftsol(:,is_sic)
     659  !GG                tice(:)=ftsol(:,is_sic)
     660                  tice_slab(:)=ftsol(:,is_sic)
    644661          ENDIF
    645662          found=phyetat0_get(seaice,"seaice","seaice",0.)
    … …  
    688705  end if
    689706
     707  !GG
     708  ! Sea ice
     709  !IF (iflag_seaice == 2) THEN
     710
     711  found=phyetat0_get(hice,"hice","Ice thickness",0.)
     712  IF (.NOT. found) THEN
     713       PRINT*, "phyetat0: Le champ <hice> est absent"
     714       PRINT*, "Initialisation a hice=1m "
     715       hice(:)=1.0
     716  END IF
     717  found=phyetat0_get(tice,"tice","Sea Ice temperature",0.)
     718  IF (.NOT. found) THEN
     719       PRINT*, "phyetat0: Le champ <tice> est absent"
     720       PRINT*, "Initialisation a tsol_sic"
     721       tice(:)=ftsol(:,is_sic)
     722  END IF
     723  found=phyetat0_get(bilg_cumul,"bilg_cumul","Flux conductivite + transmit sea-ice",0.)
     724  IF (.NOT. found) THEN
     725       PRINT*, "phyetat0: Le champ <bilg_cumul> est absent"
     726       PRINT*, "Initialisation a zero"
     727       bilg_cumul(:)=0.0
     728  END IF
     729
     730  !END IF
     731  !GG
    690732  ! on ferme le fichier
    691733  CALL close_startphy
    … …  
    694736
    695737  if ( iflag_physiq <= 1 ) then
    696   CALL pbl_surface_init(fder, snow, qsurf, tsoil)
     738  !GG CALL pbl_surface_init(fder, snow, qsurf, tsoil)
     739  CALL pbl_surface_init(fder, snow, qsurf, tsoil, hice, tice, bilg_cumul)
     740  !GG
    697741  endif
    698742

  • LMDZ6/trunk/libf/phylmd/phyredem.f90

    r5645 r5662  
    3434                                du_gwd_rando, du_gwd_front, u10m, v10m, &
    3535                                treedrg, solswfdiff, delta_sal, ds_ns, dt_ns, &
    36                                 delta_sst, ratqs_inter_, dter, dser, dt_ds,  &
    37                                 frac_tersrf, z0m_tersrf, ratio_z0m_z0h_tersrf, &
    38                                 albedo_tersrf, beta_tersrf, inertie_tersrf,  &
    39                                 hcond_tersrf, tsurfi_tersrf, tsoili_tersrf, tsoil_depth, &
    40                                 qsurf_tersrf, tsurf_tersrf, tsoil_tersrf, tsurf_new_tersrf, &
    41                                 cdragm_tersrf, cdragh_tersrf, &
    42                                 swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf
     36!GG                                delta_sst, ratqs_inter_, dter,
     37                                !dser, dt_ds
     38                                delta_sst, ratqs_inter_, dter, dser,&
     39                                & dt_ds, hice, tice, bilg_cumul, !GG
     40                                frac_tersrf, z0m_tersrf,&
     41                                     & ratio_z0m_z0h_tersrf,&
     42                                     & albedo_tersrf, beta_tersrf,&
     43                                     & inertie_tersrf, hcond_tersrf,&
     44                                     & tsurfi_tersrf, tsoili_tersrf,&
     45                                     & tsoil_depth, qsurf_tersrf,&
     46                                     & tsurf_tersrf, tsoil_tersrf,&
     47                                     & tsurf_new_tersrf,&
     48                                     & cdragm_tersrf, cdragh_tersrf,&
     49                                     & swnet_tersrf, lwnet_tersrf,&
     50                                     & fluxsens_tersrf, fluxlat_tersrf
    4351
    4452  USE geometry_mod, ONLY : longitude_deg, latitude_deg
    … …  
    4856  USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send, carbon_cycle_rad, RCO2_glo
    4957  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic, epsfra
    50   USE surface_data, ONLY: type_ocean, version_ocean
    51   USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice, fsic
     58!GG  USE surface_data, ONLY: type_ocean, version_ocean
     59  USE surface_data, ONLY: type_ocean, version_ocean, iflag_seaice, iflag_seaice_alb, &
     60                          iflag_leads
     61!GG
     62  USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice_slab, fsic
    5263  USE time_phylmdz_mod, ONLY: annee_ref, day_end, itau_phy, pdtphys
    5364  use config_ocean_skin_m, only: activate_ocean_skin
    … …  
    121132  IF (carbon_cycle_rad) tab_cntrl(17) = RCO2_glo
    122133  !PRINT*, "PC : phyredem RCO2_glo =",RCO2_glo
     134  !GG
     135  tab_cntrl(18 ) = iflag_seaice
     136  tab_cntrl(19 ) = iflag_seaice_alb
     137  tab_cntrl(20 ) = iflag_leads
     138  !GG
    123139
    124140  DO pass=1,2   ! pass=1 netcdf definition ; pass=2 netcdf write
    … …  
    225241    CALL put_field_srf1(pass,"SNOW", "Neige", snow(:,:))
    226242
     243    !GG
     244    CALL put_field(pass,"hice", "Ice thickness", hice)
     245    CALL put_field(pass,"tice", "Sea Ice temperature", tice)
     246    CALL put_field(pass,"bilg_cumul", "Flux conductivite + transmit sea-ice", bilg_cumul)
     247    !GG
     248
    227249    CALL put_field(pass,"RADS", "Rayonnement net a la surface", radsol)
    228250
    … …  
    400422        IF (version_ocean == 'sicINT') THEN
    401423            CALL put_field(pass,"seaice", "Slab seaice (kg/m2)", seaice)
    402             CALL put_field(pass,"slab_tice", "Slab sea ice temperature", tice)
     424            CALL put_field(pass,"slab_tice", "Slab sea ice temperature", tice_slab)
    403425        END IF
    404426    END IF

  • LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90

    r5649 r5662  
    176176      !$OMP THREADPRIVATE( d_xt_decroiss)
    177177#endif
     178! GG
     179! diagnostique de la glace de mer
     180      REAL, SAVE, ALLOCATABLE :: fcds(:)
     181      !$OMP THREADPRIVATE(fcds)
     182      REAL, SAVE, ALLOCATABLE :: fcdi(:)
     183      !$OMP THREADPRIVATE(fcdi)
     184      REAL, SAVE, ALLOCATABLE :: dh_basal_growth(:)
     185      !$OMP THREADPRIVATE(dh_basal_growth)
     186      REAL, SAVE, ALLOCATABLE :: dh_basal_melt(:)
     187      !$OMP THREADPRIVATE(dh_basal_melt)
     188      REAL, SAVE, ALLOCATABLE :: dh_top_melt(:)
     189      !$OMP THREADPRIVATE(dh_top_melt)
     190      REAL, SAVE, ALLOCATABLE :: dh_snow2sic(:)
     191      !$OMP THREADPRIVATE(dh_snow2sic)
     192      REAL, SAVE, ALLOCATABLE :: dtice_melt(:)
     193      !$OMP THREADPRIVATE(dtice_melt)
     194      REAL, SAVE, ALLOCATABLE :: dtice_snow2sic(:)
     195      !$OMP THREADPRIVATE(dtice_snow2sic)
     196! GG
    178197
    179198! tendance du a la conersion Ec -> E thermique
    … …  
    9971016      ALLOCATE(load_tmp9(klon))
    9981017      ALLOCATE(load_tmp10(klon))
     1018!GG
     1019      ALLOCATE(fcds(klon))
     1020      ALLOCATE(fcdi(klon))
     1021      ALLOCATE(dh_basal_growth(klon))
     1022      ALLOCATE(dh_basal_melt(klon))
     1023      ALLOCATE(dh_top_melt(klon))
     1024      ALLOCATE(dh_snow2sic(klon))
     1025      ALLOCATE(dtice_melt(klon))
     1026      ALLOCATE(dtice_snow2sic(klon))
     1027!GG
    9991028
    10001029!IM ajout variables CFMIP2/CMIP5
    … …  
    14481477      DEALLOCATE(dv_gwd_rando,dv_gwd_front)
    14491478      DEALLOCATE(east_gwstress,west_gwstress)
     1479!GG
     1480      DEALLOCATE(fcds)
     1481      DEALLOCATE(fcdi)
     1482      DEALLOCATE(dh_basal_growth)
     1483      DEALLOCATE(dh_basal_melt)
     1484      DEALLOCATE(dh_top_melt)
     1485      DEALLOCATE(dh_snow2sic)
     1486      DEALLOCATE(dtice_melt)
     1487      DEALLOCATE(dtice_snow2sic)
     1488!GG
    14501489
    14511490!IM ajout variables CFMIP2/CMIP5

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

    r5655 r5662  
    424424  TYPE(ctrl_out), SAVE :: o_fsnow = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
    425425    'fsnow', 'Surface snow area fraction', '-', (/ ('', i=1, 10) /))
     426!GG
     427  TYPE(ctrl_out), SAVE :: o_hice = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     428    'hice', 'Sea Ice Thickness', 'm', (/ ('', i=1, 10) /))
     429  TYPE(ctrl_out), SAVE :: o_fcds = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     430    'fcds', 'Cond. flux snow on sea ice', 'W/m2', (/ ('', i=1, 10) /))
     431  TYPE(ctrl_out), SAVE :: o_fcdi = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     432    'fcdi', 'Cond. flux sea ice', 'W/m2', (/ ('', i=1, 10) /))
     433  TYPE(ctrl_out), SAVE :: o_dh_basal_growth = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     434    'dh_basal_growth', 'Sea ice thickness tendency due to basal growth', 'm/day', (/ ('', i=1, 10) /))
     435  TYPE(ctrl_out), SAVE :: o_dh_basal_melt = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     436    'dh_basal_melt', 'Sea ice thickness tendency due to basal melt', 'm/day', (/ ('', i=1, 10) /))
     437  TYPE(ctrl_out), SAVE :: o_dh_top_melt = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     438    'dh_top_melt', 'Sea ice thickness tendency due to melt from above', 'm/day', (/ ('', i=1, 10) /))
     439  TYPE(ctrl_out), SAVE :: o_dh_snow2sic = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     440    'dh_snow2sic', 'Sea ice thickness tendency due snow conversion', 'm/day', (/ ('', i=1, 10) /))
     441  TYPE(ctrl_out), SAVE :: o_tice = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     442    'tice', 'Sea Ice Temperature', 'K', (/ ('', i=1, 10) /))
     443  TYPE(ctrl_out), SAVE :: o_dtice_melt = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     444    'dtice_melt', 'Sea Ice Temperature tendency due to >0 tsol', 'K/day', (/ ('', i=1, 10) /))
     445  TYPE(ctrl_out), SAVE :: o_dtice_snow2sic = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     446    'dtice_snow2sic', 'Sea Ice Temperature tendency due snow conversion', 'K/day', (/ ('', i=1, 10) /))
     447  TYPE(ctrl_out), SAVE :: o_bilg_cumul = ctrl_out((/ 1, 1, 1, 5, 10, 10, 11, 11, 11, 11/), &
     448    'bilg_cumul', 'Flux conductivite et transmis', 'W/m2', (/ ('', i=1, 10) /))
     449!GG
    426450  TYPE(ctrl_out), SAVE :: o_tops = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11, 11/), &
    427451    'tops', 'Solar rad. at TOA', 'W/m2', (/ ('', i=1, 10) /))

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

    r5655 r5662  
    7878         o_tke_srf, o_tke_max_srf,o_dltpbltke_srf, o_wstar, &
    7979         o_l_mixmin,o_l_mix, &
     80!GG
     81         o_hice, o_tice, o_bilg_cumul,&
     82!GG
    8083         o_cdrm, o_cdrh, o_cldl, o_cldm, o_cldh, &
    8184         o_cldt, o_JrNt, o_cldljn, o_cldmjn, &
    … …  
    245248         o_tks, o_taur, o_sss, &
    246249!FC
    247          o_zxfluxt,o_zxfluxq
     250!GG         o_zxfluxt,o_zxfluxq
     251         o_zxfluxt,o_zxfluxq, &
     252         o_fcds, o_fcdi, o_dh_basal_growth, o_dh_basal_melt, &
     253         o_dh_top_melt, o_dh_snow2sic, &
     254         o_dtice_melt, o_dtice_snow2sic
     255!GG
    248256
    249257#ifdef CPP_ECRAD
    … …  
    294302         v10m, pbl_tke, wake_delta_pbl_TKE, &
    295303         delta_tsurf, &
     304!GG
     305         hice, tice, bilg_cumul,&
     306!GG
    296307         wstar, cape, ema_pcb, ema_pct, &
    297308         ema_cbmf, Mipsh, Ma, fm_therm, ale_bl, alp_bl, ale, &
    … …  
    424435         zxfluxt,zxfluxq, &
    425436! offline
    426          da, mp, phi, wght_cvfd
     437!GG         da, mp, phi, wght_cvfd
     438         da, mp, phi, wght_cvfd, &
     439         fcds, fcdi, dh_basal_growth, dh_basal_melt, &
     440         dh_top_melt, dh_snow2sic, &
     441         dtice_melt, dtice_snow2sic
     442!GG
    427443    USE phys_output_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
    428444         reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra
    … …  
    467483         ok_4xCO2atm, tkt, tks, taur, sss
    468484
    469     USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice, seaice, &
     485!GG    USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice, seaice, &
     486    USE ocean_slab_mod, ONLY: nslay, tslab, slab_bilg, tice_slab, seaice, &
    470487        slab_ekman,slab_hdiff,slab_gm,dt_ekman, dt_hdiff, dt_gm, dt_qflux
    471488    USE pbl_surface_mod, ONLY: snow, ftsoil
    … …  
    475492#endif
    476493    USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
    477     USE surface_data, ONLY: type_ocean, version_ocean, ok_veget, landice_opt
     494!GG    USE surface_data, ONLY: type_ocean, version_ocean, ok_veget, landice_opt
     495    USE surface_data, ONLY: type_ocean, version_ocean, ok_veget, landice_opt, &
     496        iflag_seaice, iflag_seaice_alb
     497    !GG
    478498    USE aero_mod, ONLY: naero_tot, id_STRAT_phy
    479499    USE ioipsl, ONLY: histend, histsync
    … …  
    17391759          CALL histwrite_phy(o_alp_bl_stat, alp_bl_stat)
    17401760       ENDIF  !(iflag_clos_bl>=1)
     1761!GG
     1762       IF (iflag_seaice>0 ) THEN
     1763       !IF (ok_hice ) THEN
     1764          CALL histwrite_phy(o_hice, hice)
     1765          CALL histwrite_phy(o_tice, tice)
     1766          CALL histwrite_phy(o_bilg_cumul, bilg_cumul)
     1767          CALL histwrite_phy(o_fcds, fcds)
     1768          CALL histwrite_phy(o_fcdi, fcdi)
     1769          CALL histwrite_phy(o_dh_basal_growth, dh_basal_growth)
     1770          CALL histwrite_phy(o_dh_basal_melt, dh_basal_melt)
     1771          CALL histwrite_phy(o_dh_top_melt, dh_top_melt)
     1772          CALL histwrite_phy(o_dh_snow2sic, dh_snow2sic)
     1773          CALL histwrite_phy(o_dtice_melt, dtice_melt)
     1774          CALL histwrite_phy(o_dtice_snow2sic, dtice_snow2sic)
     1775       END IF
     1776!GG
    17411777!!! fin nrlmd le 10/04/2012
    17421778       ! Output of slab ocean variables
    … …  
    17541790          IF (version_ocean=='sicINT') THEN
    17551791              CALL histwrite_phy(o_slab_bilg, slab_bilg)
    1756               CALL histwrite_phy(o_slab_tice, tice)
     1792              CALL histwrite_phy(o_slab_tice, tice_slab)
     1793      !GG        CALL histwrite_phy(o_slab_tice, tice)
    17571794              CALL histwrite_phy(o_slab_sic, seaice)
    17581795          ENDIF

  • LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90

    r5627 r5662  
    218218      REAL,ALLOCATABLE,SAVE :: O3sumSTD(:,:,:), O3daysumSTD(:,:,:)
    219219!$OMP THREADPRIVATE(O3sumSTD,O3daysumSTD)
     220!GG
     221      REAL,ALLOCATABLE,SAVE :: hice(:)
     222!$OMP THREADPRIVATE(hice)
     223      REAL,ALLOCATABLE,SAVE :: tice(:)
     224!$OMP THREADPRIVATE(tice)
     225      REAL,ALLOCATABLE,SAVE :: bilg_cumul(:)
     226!$OMP THREADPRIVATE(bilg_cumul)
     227!GG
    220228!IM begin
    221229      REAL,ALLOCATABLE,SAVE :: wlevSTD(:,:), ulevSTD(:,:), vlevSTD(:,:)
    … …  
    719727      ALLOCATE(zuthe(klon),zvthe(klon))
    720728      ALLOCATE(alb_neig(klon))
     729!GG
     730      ALLOCATE(hice(klon))
     731      ALLOCATE(tice(klon))
     732      ALLOCATE(bilg_cumul(klon))
     733!GG
    721734!cloud base mass flux
    722735      ALLOCATE(ema_cbmf(klon))
    … …  
    945958      DEALLOCATE(zuthe, zvthe)
    946959      DEALLOCATE(alb_neig)
     960!GG
     961      DEALLOCATE(hice)
     962      DEALLOCATE(tice)
     963      DEALLOCATE(bilg_cumul)
     964!GG
    947965      DEALLOCATE(ema_cbmf)
    948966      DEALLOCATE(ema_pcb, ema_pct)

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r5659 r5662  
    361361       rneb,  &
    362362       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
    363        zxfluxt,zxfluxq
     363!GG       zxfluxt,zxfluxq
     364       zxfluxt,zxfluxq, &
     365       fcds, fcdi, dh_basal_growth, dh_basal_melt, &
     366       dh_top_melt, dh_snow2sic, &
     367       dtice_melt, dtice_snow2sic
     368!GG
    364369       !
    365370       USE phys_local_var_mod, ONLY: zfice, dNovrN, ptconv
    … …  
    28792884            debut,     lafin, &
    28802885            longitude_deg, latitude_deg, rugoro,  zrmu0,      &
    2881             sollwdown,    cldt,      &
     2886        !GG    sollwdown,    cldt,      &
     2887            sollwdown, pphi,   cldt,      &
     2888        !GG
    28822889            rain_fall, snow_fall, bs_fall, solsw,   solswfdiff, sollw,     &
    28832890            gustiness,                                &
    … …  
    29242931            wake_delta_pbl_TKE, &
    29252932                                !>nrlmd+jyg
    2926              treedrg, &
     2933!GG             treedrg )
     2934            treedrg,hice, tice, bilg_cumul, &
     2935            fcds, fcdi, dh_basal_growth, dh_basal_melt, &
     2936            dh_top_melt, dh_snow2sic, &
     2937            dtice_melt, dtice_snow2sic , &
     2938!GG
    29272939!FC
    29282940!AM

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

    r5285 r5662  
    2121       tsurf_new, dflux_s, dflux_l, lmt_bils, &
    2222       flux_u1, flux_v1, delta_sst, delta_sal, ds_ns, dt_ns, dter, dser, &
    23        dt_ds, tkt, tks, taur, sss &
     23!GG       dt_ds, tkt, tks, taur, sss)
     24       dt_ds, tkt, tks, taur, sss, &
     25       dthetadz300,Ampl      &
     26!GG
    2427#ifdef ISO
    2528        &       ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
    … …  
    127130    ! (tks / tkt) * dTer, in K
    128131
     132!GG
     133    REAL, DIMENSION(klon), INTENT(IN)        :: dthetadz300
     134    REAL, DIMENSION(klon), INTENT(OUT)        :: Ampl
     135!
     136
    129137    ! Output variables
    130138    !**************************************************************************
    … …  
    270278            radsol, snow, agesno, &
    271279            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    272             tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa &
     280!GG           tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa)
     281            tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa, &
     282            dthetadz300,  pctsrf, Ampl &
     283!GG
    273284#ifdef ISO
    274285            ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &

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

    r5285 r5662  
    2121       z0m, z0h, SFRWL, alb_dir_new, alb_dif_new, evap, fluxsens, fluxlat, & 
    2222       tsurf_new, dflux_s, dflux_l, &
    23        flux_u1, flux_v1 &
     23!GG       flux_u1, flux_v1)
     24       flux_u1, flux_v1, hice, tice,bilg_cumul, &
     25       fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &
     26       dtice_melt, dtice_snow2sic &
     27!GG
    2428#ifdef ISO
    2529         &      ,xtprecip_rain, xtprecip_snow,xtspechum,Roce, &
    … …  
    101105    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l     
    102106    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
     107!GG
     108    REAL, DIMENSION(klon), INTENT(INOUT)       :: hice, tice, bilg_cumul
     109    REAL, DIMENSION(klon), INTENT(INOUT)       :: fcds,fcdi, dh_basal_growth,dh_basal_melt
     110    REAL, DIMENSION(klon), INTENT(INOUT)       :: dh_top_melt, dh_snow2sic, dtice_melt, dtice_snow2sic
     111!GG
    103112#ifdef ISO
    104113    REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap
    … …  
    169178            AcoefH, AcoefQ, BcoefH, BcoefQ, &
    170179            AcoefU, AcoefV, BcoefU, BcoefV, &
    171             ps, u1, v1, gustiness, &
     180!GG            ps, u1, v1, gustiness, &
     181            ps, u1, v1, gustiness,pctsrf, &
     182!GG
    172183            radsol, snow, qsol, agesno, tsoil, &
    173184            qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    174             tsurf_new, dflux_s, dflux_l, rhoa &
     185!GG            tsurf_new, dflux_s, dflux_l, rhoa)
     186            tsurf_new, dflux_s, dflux_l,rhoa,swnet,hice, tice, bilg_cumul, &
     187            fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, &
     188            dtice_melt, dtice_snow2sic &
     189!GG
    175190#ifdef ISO
    176191            ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &

  • LMDZ6/trunk/libf/phylmd/surface_data.f90

    r5268 r5662  
    1717  CHARACTER(len=6), SAVE :: type_ocean    ! force/slab/couple
    1818  !$OMP THREADPRIVATE(type_ocean)
     19
     20  !GG
     21  INTEGER, SAVE          :: iflag_seaice
     22  !$OMP THREADPRIVATE(iflag_seaice)
     23  INTEGER, SAVE          :: iflag_seaice_alb
     24  !$OMP THREADPRIVATE(iflag_seaice_alb)
     25  INTEGER, SAVE          :: iflag_leads
     26  !$OMP THREADPRIVATE(iflag_leads)
     27
     28  !For sea-ice
     29  REAL,SAVE             :: sice_cond
     30  !$OMP THREADPRIVATE(sice_cond)
     31  REAL,SAVE             :: sisno_cond
     32  !$OMP THREADPRIVATE(sisno_cond)
     33  REAL,SAVE             :: sisno_den
     34  !$OMP THREADPRIVATE(sisno_den)
     35  REAL,SAVE             :: sisno_min
     36  !$OMP THREADPRIVATE(sisno_min)
     37  REAL,SAVE             :: sithick_min
     38  !$OMP THREADPRIVATE(sithick_min)
     39  REAL,SAVE             :: sisno_wfact
     40  !$OMP THREADPRIVATE(sisno_wfact)
     41  REAL,SAVE             :: amax_n
     42  !$OMP THREADPRIVATE(amax_n)
     43  REAL,SAVE             :: amax_s
     44  !$OMP THREADPRIVATE(amax_s)
     45  REAL,SAVE             :: rn_alb_sdry
     46  !$OMP THREADPRIVATE(rn_alb_sdry)
     47  REAL,SAVE             :: rn_alb_smlt
     48  !$OMP THREADPRIVATE(rn_alb_smlt)
     49  REAL,SAVE             :: rn_alb_idry
     50  !$OMP THREADPRIVATE(rn_alb_idry)
     51  REAL,SAVE             :: rn_alb_imlt
     52  !$OMP THREADPRIVATE(rn_alb_imlt)
     53  REAL,SAVE             :: si_pen_frac
     54  !$OMP THREADPRIVATE(si_pen_frac)
     55  REAL,SAVE             :: si_pen_ext
     56  !$OMP THREADPRIVATE(si_pen_ext)
     57  REAL,SAVE             :: fseaN
     58  !$OMP THREADPRIVATE(fseaN)
     59  REAL,SAVE             :: fseaS
     60  !$OMP THREADPRIVATE(fseaS)
     61  !GG
    1962
    2063  ! if type_ocean=couple : version_ocean=opa8 ou nemo

  • LMDZ6/trunk/libf/phylmdiso/limit_read_mod.F90

    r5217 r5662  
    2222  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: sst
    2323!$OMP THREADPRIVATE(sst)   
     24!GG
     25  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: sih
     26!$OMP THREADPRIVATE(sih)
     27!GG
    2428#ifdef ISO
    2529  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: tuoce
    … …  
    254258  END SUBROUTINE limit_read_sst
    255259
     260!GG
     261  SUBROUTINE limit_read_hice(knon, knindex, hice_out)
     262!
     263! This subroutine returns the sea surface temperature already read from limit.nc.
     264!
     265    USE dimphy, ONLY : klon
     266
     267    INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
     268    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
     269    REAL, DIMENSION(klon), INTENT(OUT)   :: hice_out
     270
     271    INTEGER :: i
     272
     273    DO i = 1, knon
     274       hice_out(i) = sih(knindex(i))
     275    END DO
     276
     277  END SUBROUTINE limit_read_hice
     278!GG
    256279!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    257280!!
    … …  
    273296    USE mod_grid_phy_lmdz
    274297    USE mod_phys_lmdz_para
    275     USE surface_data, ONLY : type_ocean, ok_veget
     298    !GG USE surface_data, ONLY : type_ocean, ok_veget
     299    USE surface_data, ONLY : type_ocean, ok_veget, iflag_seaice, amax_n, amax_s
     300    !GG
    276301    USE netcdf
    277302    USE indice_sol_mod
    … …  
    284309    USE phys_cal_mod, ONLY : calend, year_len
    285310    USE print_control_mod, ONLY: lunout, prt_level
    286     USE lmdz_XIOS, ONLY: xios_recv_field
     311    USE lmdz_xios, ONLY: xios_recv_field, using_xios
    287312   
    288313    IMPLICIT NONE
    … …  
    319344    REAL, DIMENSION(klon_mpi)                 :: rug_mpi  ! rugosity at global grid
    320345    REAL, DIMENSION(klon_mpi)                 :: alb_mpi  ! albedo at global grid
     346!GG
     347    REAL, DIMENSION(klon_glo)                 :: sih_glo  ! albedo at global grid
     348    REAL, DIMENSION(klon_mpi)                 :: sih_mpi  ! albedo at global grid
     349!GG
    321350
    322351    CHARACTER(len=20)                         :: modname='limit_read_mod'     
    … …  
    344373       END IF
    345374
     375       !GG
     376       IF (iflag_seaice==1) THEN
     377             ALLOCATE(sih(klon), stat=ierr)
     378             IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating sih',1)
     379       ENDIF
     380       !GG
     381
    346382       IF ( .NOT. ok_veget ) THEN
    347383          ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
    … …  
    422458         IF ( type_ocean /= 'couple') THEN                   
    423459             IF (is_omp_master) CALL xios_recv_field("sst_limin",sst_mpi)
     460             !GG
     461             IF (is_omp_master) CALL xios_recv_field("sih_limin",sih_mpi)
     462             !GG
    424463         ENDIF
    425464       
    … …  
    431470       IF ( type_ocean /= 'couple') THEN
    432471          CALL Scatter_omp(sst_mpi,sst)
     472          !GG
     473          CALL Scatter_omp(sih_mpi,sih)
     474          !GG
    433475          CALL Scatter_omp(pct_mpi(:,is_oce),pctsrf(:,is_oce))
    434476          CALL Scatter_omp(pct_mpi(:,is_sic),pctsrf(:,is_sic))
    … …  
    481523             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
    482524
     525! GG
     526! Account for leads
     527             IF (iflag_seaice>0) THEN
     528               DO ii=1,klon_glo/2
     529                 if (pct_glo(ii,is_sic)>amax_n) THEN
     530                    pct_glo(ii,is_oce)=pct_glo(ii,is_oce)+(pct_glo(ii,is_sic)-amax_n)
     531                    pct_glo(ii,is_sic)=amax_n
     532                 end if
     533               ENDDO
     534               DO ii=klon_glo/2,klon_glo
     535               if (pct_glo(ii,is_sic)>amax_s) THEN
     536                    pct_glo(ii,is_oce)=pct_glo(ii,is_oce)+(pct_glo(ii,is_sic)-amax_s)
     537                    pct_glo(ii,is_sic)=amax_s
     538               end if
     539               ENDDO
     540             ENDIF
     541!GG
    483542
    484543! Read land and continentals fraction only if asked for
    … …  
    513572             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
    514573             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
    515          
     574!GG
     575             IF (iflag_seaice == 1) THEN
     576               ierr = NF90_INQ_VARID(nid, 'HICE', nvarid)
     577               IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <HICE> est absent',1)
     578
     579               ierr = NF90_GET_VAR(nid,nvarid,sih_glo(:),start,epais)
     580               IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <HICE>' ,1)
     581             ENDIF
     582            !GG
     583         
    516584#ifdef ISO
    517585             IF ((iso_HTO.gt.0).and.(ok_prod_nucl_tritium)) THEN
    … …  
    572640       IF ( type_ocean /= 'couple') THEN
    573641          CALL Scatter(sst_glo,sst)
     642          !GG
     643          IF (iflag_seaice==1) THEN
     644             CALL Scatter(sih_glo,sih)
     645          END IF
     646          !GG
    574647          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
    575648          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))

  • LMDZ6/trunk/libf/phylmdiso/phyaqua_mod.F90

    r5645 r5662  
    351351
    352352
    353     CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
     353!GG
     354    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil, hice, tice, bilg_cumul)
     355!    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
     356!GG
    354357#ifdef ISO
    355358  CALL pbl_surface_init_iso(xtsnsrf,Rland_ice)

  • LMDZ6/trunk/libf/phylmdiso/phyetat0_mod.F90

    r5645 r5662  
    2020#endif
    2121  USE phyetat0_get_mod, ONLY : phyetat0_get, phyetat0_srf
    22   USE surface_data,     ONLY : type_ocean, version_ocean
     22!GG  USE surface_data,     ONLY : type_ocean, version_ocean
     23  USE surface_data,     ONLY : type_ocean, version_ocean, iflag_seaice, &
     24                                   iflag_seaice_alb, iflag_leads
     25!GG
    2326  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, phys_tstep, &
    2427       qsol, fevap, z0m, z0h, agesno, &
    … …  
    3740       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
    3841       ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, dter, dser, &
    39        dt_ds, ratqs_inter_
     42!GG       dt_ds, ratqs_inter_
     43       dt_ds, ratqs_inter_, &
     44       hice, tice, bilg_cumul
     45!GG
     46
    4047!FC
    4148  USE geometry_mod,     ONLY: longitude_deg, latitude_deg
    … …  
    4653  USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl, carbon_cycle_tr, carbon_cycle_rad, co2_send, RCO2_glo
    4754  USE indice_sol_mod,   ONLY: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
    48   USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice, ocean_slab_init
     55  !GG USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice, ocean_slab_init
     56  USE ocean_slab_mod,   ONLY: nslay, tslab, seaice, tice_slab, ocean_slab_init
     57  !GG
    4958  USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
    5059  use wxios_mod, ONLY: missing_val_xios => missing_val, using_xios
    … …  
    179188  IF (ok_orolf) tab_cntrl(11) =1.
    180189  IF (ok_limitvrai) tab_cntrl(12) =1.
     190  !GG
     191  tab_cntrl(18) =iflag_seaice
     192  tab_cntrl(19) =iflag_seaice_alb
     193  tab_cntrl(20) =iflag_leads
     194  !GG
    181195
    182196  itau_phy = tab_cntrl(15)
    … …  
    624638      ! Sea ice variables
    625639      IF (version_ocean == 'sicINT') THEN
    626           found=phyetat0_get(tice,"slab_tice","slab_tice",0.)
     640          found=phyetat0_get(tice_slab,"slab_tice","slab_tice",0.)
     641  !GG        found=phyetat0_get(tice,"slab_tice","slab_tice",0.)
    627642          IF (.NOT. found) THEN
    628               PRINT*, "phyetat0: Le champ <tice> est absent"
     643  !GG            PRINT*, "phyetat0: Le champ <tice> est absent"
     644              PRINT*, "phyetat0: Le champ <tice_slab> est absent"
    629645              PRINT*, "Initialisation a tsol_sic"
    630                   tice(:)=ftsol(:,is_sic)
     646  !GG                tice(:)=ftsol(:,is_sic)
     647                  tice_slab(:)=ftsol(:,is_sic)
    631648          ENDIF
    632649          found=phyetat0_get(seaice,"seaice","seaice",0.)
    … …  
    675692  end if
    676693
     694  !GG
     695  ! Sea ice
     696  !IF (iflag_seaice == 2) THEN
     697
     698  found=phyetat0_get(hice,"hice","Ice thickness",0.)
     699  IF (.NOT. found) THEN
     700       PRINT*, "phyetat0: Le champ <hice> est absent"
     701       PRINT*, "Initialisation a hice=1m "
     702       hice(:)=1.0
     703  END IF
     704  found=phyetat0_get(tice,"tice","Sea Ice temperature",0.)
     705  IF (.NOT. found) THEN
     706       PRINT*, "phyetat0: Le champ <tice> est absent"
     707       PRINT*, "Initialisation a tsol_sic"
     708       tice(:)=ftsol(:,is_sic)
     709  END IF
     710  found=phyetat0_get(bilg_cumul,"bilg_cumul","Flux conductivite + transmit sea-ice",0.)
     711  IF (.NOT. found) THEN
     712       PRINT*, "phyetat0: Le champ <bilg_cumul> est absent"
     713       PRINT*, "Initialisation a zero"
     714       bilg_cumul(:)=0.0
     715  END IF
     716
     717  !END IF
     718  !GG
    677719  ! on ferme le fichier
    678720  CALL close_startphy
    … …  
    686728!#endif
    687729  if ( iflag_physiq <= 1 ) then
    688   CALL pbl_surface_init(fder, snow, qsurf, tsoil)
     730  !GG CALL pbl_surface_init(fder, snow, qsurf, tsoil)
     731  CALL pbl_surface_init(fder, snow, qsurf, tsoil, hice, tice, bilg_cumul)
     732  !GG
    689733#ifdef ISO
    690734  CALL pbl_surface_init_iso(xtsnow,Rland_ice)

  • LMDZ6/trunk/libf/phylmdiso/phyredem.F90

    r5645 r5662  
    3131                                du_gwd_rando, du_gwd_front, u10m, v10m, &
    3232                                treedrg, solswfdiff, delta_sal, ds_ns, dt_ns, &
    33                                 delta_sst, ratqs_inter_, dter, dser, dt_ds
     33!GG                                delta_sst, ratqs_inter_, dter, dser, dt_ds
     34                                delta_sst, ratqs_inter_, dter, dser, dt_ds, &
     35                                hice, tice, bilg_cumul
     36!GG
    3437#ifdef ISO
    3538  USE phys_state_var_mod, ONLY: xtsol, fxtevap,xtrain_fall, xtsnow_fall,     &
    … …  
    5154  USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send, carbon_cycle_rad, RCO2_glo
    5255  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic, epsfra
    53   USE surface_data, ONLY: type_ocean, version_ocean
    54   USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice, fsic
     56!GG  USE surface_data, ONLY: type_ocean, version_ocean
     57  USE surface_data, ONLY: type_ocean, version_ocean, iflag_seaice, iflag_seaice_alb, &
     58                          iflag_leads
     59!GG
     60  USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice_slab, fsic
    5561  USE time_phylmdz_mod, ONLY: annee_ref, day_end, itau_phy, pdtphys
    5662  use config_ocean_skin_m, only: activate_ocean_skin
    … …  
    136142  IF (carbon_cycle_rad) tab_cntrl(17) = RCO2_glo
    137143  !PRINT*, "PC : phyredem RCO2_glo =",RCO2_glo
     144  !GG
     145  tab_cntrl(18 ) = iflag_seaice
     146  tab_cntrl(19 ) = iflag_seaice_alb
     147  tab_cntrl(20 ) = iflag_leads
     148  !GG
    138149
    139150  DO pass=1,2   ! pass=1 netcdf definition ; pass=2 netcdf write
    … …  
    214225
    215226    CALL put_field_srf1(pass,"SNOW", "Neige", snow(:,:))
     227
     228    !GG
     229    CALL put_field(pass,"hice", "Ice thickness", hice)
     230    CALL put_field(pass,"tice", "Sea Ice temperature", tice)
     231    CALL put_field(pass,"bilg_cumul", "Flux conductivite + transmit sea-ice", bilg_cumul)
     232    !GG
    216233
    217234    CALL put_field(pass,"RADS", "Rayonnement net a la surface", radsol)
    … …  
    388405        IF (version_ocean == 'sicINT') THEN
    389406            CALL put_field(pass,"seaice", "Slab seaice (kg/m2)", seaice)
    390             CALL put_field(pass,"slab_tice", "Slab sea ice temperature", tice)
     407            CALL put_field(pass,"slab_tice", "Slab sea ice temperature", tice_slab)
    391408        END IF
    392409    END IF

  • LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

    r5659 r5662  
    402402       rneb,  &
    403403       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
    404        zxfluxt,zxfluxq
     404!GG       zxfluxt,zxfluxq
     405       zxfluxt,zxfluxq, &
     406       fcds, fcdi, dh_basal_growth, dh_basal_melt, &
     407       dh_top_melt, dh_snow2sic, &
     408       dtice_melt, dtice_snow2sic
     409!GG
    405410
    406411
    … …  
    32993304            debut,     lafin, &
    33003305            longitude_deg, latitude_deg, rugoro,  zrmu0,      &
    3301             sollwdown,    cldt,      &
     3306        !GG    sollwdown,    cldt,      &
     3307            sollwdown, pphi,   cldt,      &
     3308        !GG
    33023309            rain_fall, snow_fall, bs_fall, solsw,   solswfdiff, sollw,     &
    33033310            gustiness,                                &
    … …  
    33473354            wake_delta_pbl_TKE, &
    33483355                                !>nrlmd+jyg
    3349              treedrg &           
     3356!GG             treedrg )
     3357            treedrg,hice, tice, bilg_cumul, &
     3358            fcds, fcdi, dh_basal_growth, dh_basal_melt, &
     3359            dh_top_melt, dh_snow2sic, &
     3360            dtice_melt, dtice_snow2sic , &
     3361!GG
    33503362!AM
    3351             , tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
     3363            tsurf_tersrf, tsoil_tersrf, qsurf_tersrf, tsurf_new_tersrf, &
    33523364            cdragm_tersrf, cdragh_tersrf, &
    33533365            swnet_tersrf, lwnet_tersrf, fluxsens_tersrf, fluxlat_tersrf &
Note: See TracChangeset for help on using the changeset viewer.