New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6333 – NEMO

Changeset 6333


Ignore:
Timestamp:
2016-02-19T13:11:22+01:00 (8 years ago)
Author:
timgraham
Message:

Merged in revision r6324 of v3_6_STABLE. Also added in changes from change set at r5540 as these were somehow missed

Location:
branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM
Files:
1 deleted
75 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/ARCH/CMCC/arch-ifort_athena_xios.fcm

    r6237 r6333  
    3434 
    3535# required modules 
    36 # module load  INTEL/intel_xe_2013 NETCDF/netcdf-4.3_parallel NETCDF/parallel-netcdf-1.3.1 HDF5/hdf5-1.8.11_parallel 
     36# module load  INTEL/intel_xe_2013 NETCDF/netcdf-4.3_parallel_shared NETCDF/parallel-netcdf-1.3.1 HDF5/hdf5-1.8.11_parallel_shared 
    3737 
    3838# Environment variables set by user. Others should automatically define when loading modules. 
    39 # export XIOS=/users/home/models/nemo/xios 
     39#export XIOS=/users/home/models/nemo/xios 
     40#export HDF5=/users/home/opt/hdf5/hdf5-1.8.11_parallel_shared 
     41#export NETCDF=/users/home/opt/netcdf/netcdf-4.3_parallel_shared 
    4042 
    41 %NCDF_INC            -I$NETCDF/include -I$PNETCDF/include 
    42 %NCDF_LIB            -L$NETCDF/lib -lnetcdff -lnetcdf -L$PNETCDF/lib -lpnetcdf 
    43 %HDF5_INC            -I$PHDF5/include 
    44 %HDF5_LIB            -L$PHDF5/lib -lhdf5_hl -lhdf5 
    45 %XIOS_INC            -I$XIOS/inc 
    46 %XIOS_LIB            -L$XIOS/lib -lxios 
     43%NCDF_INC            -I${NETCDF}/include  
     44%NCDF_LIB            -L${NETCDF}/lib -lnetcdff -lnetcdf 
     45%HDF5_INC            -I${HDF5}/include 
     46%HDF5_LIB            -L${HDF5}/lib -lhdf5_hl -lhdf5 
     47%XIOS_INC            -I${XIOS}/inc 
     48%XIOS_LIB            -L${XIOS}/lib -lxios 
    4749%CPP                 cpp 
    4850%FC                  mpiifort 
    49 %FCFLAGS             -r8 -O3 -xHost -fp-model source -traceback 
     51%FCFLAGS             -r8 -O3 -xHost -fp-model source -traceback ${CFLAGS}  
    5052%FFLAGS              %FCFLAGS 
    5153%LD                  mpiifort 
    5254%FPPFLAGS            -P -C -traditional 
    53 %LDFLAGS             -lstdc++ -lz -lgpfs -lcurl 
     55%LDFLAGS             -lstdc++ -lz -lgpfs -lcurl  ${LDFLAGS} 
    5456%AR                  ar  
    5557%ARFLAGS             -r 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg

    r5501 r6333  
    200200/ 
    201201!----------------------------------------------------------------------- 
    202 &namobc        !   open boundaries parameters                           ("key_obc") 
    203 !----------------------------------------------------------------------- 
    204 / 
    205 !----------------------------------------------------------------------- 
    206202&namagrif      !  AGRIF zoom                                            ("key_agrif") 
    207203!----------------------------------------------------------------------- 
     
    369365/ 
    370366!----------------------------------------------------------------------- 
     367&namzdf_tmx_new !  new tidal mixing parameterization                    ("key_zdftmx_new") 
     368!----------------------------------------------------------------------- 
     369/ 
     370!----------------------------------------------------------------------- 
    371371&namsol        !   elliptic solver / island / free surface 
    372372!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg

    r6237 r6333  
    179179/ 
    180180!----------------------------------------------------------------------- 
    181 &namobc        !   open boundaries parameters                           ("key_obc") 
    182 !----------------------------------------------------------------------- 
    183 / 
    184 !----------------------------------------------------------------------- 
    185181&namagrif      !  AGRIF zoom                                            ("key_agrif") 
    186182!----------------------------------------------------------------------- 
     
    307303/ 
    308304!----------------------------------------------------------------------- 
     305&namzdf_tmx_new !  new tidal mixing parameterization                    ("key_zdftmx_new") 
     306!----------------------------------------------------------------------- 
     307/ 
     308!----------------------------------------------------------------------- 
    309309&namsol        !   elliptic solver / island / free surface 
    310310!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/GYRE/EXP00/namelist_cfg

    r5407 r6333  
    160160/ 
    161161!----------------------------------------------------------------------- 
    162 &namobc        !   open boundaries parameters                           ("key_obc") 
    163 !----------------------------------------------------------------------- 
    164 / 
    165 !----------------------------------------------------------------------- 
    166162&namagrif      !  AGRIF zoom                                            ("key_agrif") 
    167163!----------------------------------------------------------------------- 
     
    304300/ 
    305301!----------------------------------------------------------------------- 
     302&namzdf_tmx_new !  new tidal mixing parameterization                    ("key_zdftmx_new") 
     303!----------------------------------------------------------------------- 
     304/ 
     305!----------------------------------------------------------------------- 
    306306&namsol        !   elliptic solver / island / free surface 
    307307!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg

    r5407 r6333  
    165165/ 
    166166!----------------------------------------------------------------------- 
    167 &namobc        !   open boundaries parameters                           ("key_obc") 
    168 !----------------------------------------------------------------------- 
    169 / 
    170 !----------------------------------------------------------------------- 
    171167&namagrif      !  AGRIF zoom                                            ("key_agrif") 
    172168!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/GYRE_XIOS/EXP00/namelist_cfg

    r5407 r6333  
    154154/ 
    155155!----------------------------------------------------------------------- 
    156 &namobc        !   open boundaries parameters                           ("key_obc") 
    157 !----------------------------------------------------------------------- 
    158 / 
    159 !----------------------------------------------------------------------- 
    160156&namagrif      !  AGRIF zoom                                            ("key_agrif") 
    161157!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg

    r4990 r6333  
    165165/ 
    166166!----------------------------------------------------------------------- 
     167&namzdf_tmx_new !  new tidal mixing parameterization                    ("key_zdftmx_new") 
     168!----------------------------------------------------------------------- 
     169/ 
     170!----------------------------------------------------------------------- 
    167171&namsol        !   elliptic solver / island / free surface 
    168172!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml

    r5517 r6333  
    6161   <file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc --> 
    6262     <field field_ref="empmr"        name="wfo"      /> 
     63          <field field_ref="emp_oce"      name="emp_oce"  long_name="Evap minus Precip over ocean"              /> 
     64          <field field_ref="emp_ice"      name="emp_ice"  long_name="Evap minus Precip over ice"                /> 
    6365     <field field_ref="qsr_oce"      name="qsr_oce"  /> 
    6466     <field field_ref="qns_oce"      name="qns_oce"  /> 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_cfg

    r4995 r6333  
    168168/ 
    169169!----------------------------------------------------------------------- 
     170&namzdf_tmx_new !  new tidal mixing parameterization                    ("key_zdftmx_new") 
     171!----------------------------------------------------------------------- 
     172/ 
     173!----------------------------------------------------------------------- 
    170174&namsol        !   elliptic solver / island / free surface 
    171175!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/ORCA2_LIM_CFC_C14b/EXP00/1_namelist_cfg

    r5407 r6333  
    55!!                                    namsbc_cpl, namtra_qsr, namsbc_rnf,  
    66!!                                    namsbc_apr, namsbc_ssr, namsbc_alb) 
    7 !!              4 - lateral boundary (namlbc, namcla, namobc, namagrif, nambdy, nambdy_tide) 
     7!!              4 - lateral boundary (namlbc, namcla, namagrif, nambdy, nambdy_tide) 
    88!!              5 - bottom  boundary (nambfr, nambbc, nambbl) 
    99!!              6 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_dmp) 
     
    303303!!   namlbc        lateral momentum boundary condition 
    304304!!   namcla        cross land advection 
    305 !!   namobc        open boundaries parameters                           ("key_obc") 
    306305!!   namagrif      agrif nested grid ( read by child model only )       ("key_agrif")  
    307306!!   nambdy        Unstructured open boundaries                         ("key_bdy") 
     
    319318!----------------------------------------------------------------------- 
    320319   nn_cla      =    0      !  advection between 2 ocean pts separates by land 
    321 / 
    322 !----------------------------------------------------------------------- 
    323 &namobc        !   open boundaries parameters                           ("key_obc") 
    324 !----------------------------------------------------------------------- 
    325    ln_obc_clim = .false.   !  climatological obc data files (T) or not (F) 
    326    ln_vol_cst  = .true.    !  impose the total volume conservation (T) or not (F) 
    327    ln_obc_fla  = .false.   !  Flather open boundary condition  
    328    nn_obcdta   =    1      !  = 0 the obc data are equal to the initial state 
    329                            !  = 1 the obc data are read in 'obc.dta' files 
    330    cn_obcdta   = 'annual'  !  set to annual if obc datafile hold 1 year of data 
    331                            !  set to monthly if obc datafile hold 1 month of data 
    332    rn_dpein    =    1.     !  damping time scale for inflow at east  open boundary 
    333    rn_dpwin    =    1.     !     -           -         -       west    -      - 
    334    rn_dpnin    =    1.     !     -           -         -       north   -      - 
    335    rn_dpsin    =    1.     !     -           -         -       south   -      - 
    336    rn_dpeob    = 3000.     !  time relaxation (days) for the east  open boundary 
    337    rn_dpwob    =   15.     !     -           -         -     west    -      - 
    338    rn_dpnob    = 3000.     !     -           -         -     north   -      - 
    339    rn_dpsob    =   15.     !     -           -         -     south   -      - 
    340    rn_volemp   =    1.     !  = 0 the total volume change with the surface flux (E-P-R) 
    341                            !  = 1 the total volume remains constant 
    342320/ 
    343321!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/ORCA2_LIM_CFC_C14b/EXP00/namelist_cfg

    r5407 r6333  
    136136/ 
    137137!----------------------------------------------------------------------- 
    138 &namobc        !   open boundaries parameters                           ("key_obc") 
    139 !----------------------------------------------------------------------- 
    140 / 
    141 !----------------------------------------------------------------------- 
    142138&namagrif      !  AGRIF zoom                                            ("key_agrif") 
    143139!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist_cfg

    r4370 r6333  
    165165/ 
    166166!----------------------------------------------------------------------- 
     167&namzdf_tmx_new !  new tidal mixing parameterization                    ("key_zdftmx_new") 
     168!----------------------------------------------------------------------- 
     169/ 
     170!----------------------------------------------------------------------- 
    167171&namsol        !   elliptic solver / island / free surface 
    168172!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/SHARED/field_def.xml

    r5517 r6333  
    2323      <field_group id="grid_T" grid_ref="grid_T_2D" > 
    2424         <field id="e3t"          long_name="T-cell thickness"   standard_name="cell_thickness"   unit="m"   grid_ref="grid_T_3D"/> 
     25         <field id="e3t_0"        long_name="Initial T-cell thickness"   standard_name="ref_cell_thickness"   unit="m"   grid_ref="grid_T_3D"/> 
    2526 
    2627         <field id="toce"         long_name="temperature"         standard_name="sea_water_potential_temperature"   unit="degC"     grid_ref="grid_T_3D"/> 
     
    5960         <field id="alpha"        long_name="thermal expansion"                                                         unit="degC-1" grid_ref="grid_T_3D" /> 
    6061         <field id="beta"         long_name="haline contraction"                                                        unit="1e3"    grid_ref="grid_T_3D" /> 
    61          <field id="bn2"          long_name="squared Brunt-Vaisala frequency"                                           unit="s-1"    grid_ref="grid_T_3D" /> 
    6262         <field id="rhop"         long_name="potential density (sigma0)"        standard_name="sea_water_sigma_theta"   unit="kg/m3"  grid_ref="grid_T_3D" /> 
    6363 
     
    174174      <field_group id="SBC" grid_ref="grid_T_2D" > <!-- time step automaticaly defined based on nn_fsbc --> 
    175175         <field id="empmr"        long_name="Net Upward Water Flux"                standard_name="water_flux_out_of_sea_ice_and_sea_water"                              unit="kg/m2/s"   /> 
     176         <field id="emp_oce"      long_name="Evap minus Precip over ocean"         standard_name="evap_minus_precip_over_sea_water"                                     unit="kg/m2/s"   /> 
     177         <field id="emp_ice"      long_name="Evap minus Precip over ice"           standard_name="evap_minus_precip_over_sea_ice"                                       unit="kg/m2/s"   /> 
    176178         <field id="saltflx"      long_name="Downward salt flux"                                                                                                        unit="1e-3/m2/s" /> 
    177179         <field id="fmmflx"       long_name="Water flux due to freezing/melting"                                                                                        unit="kg/m2/s"   /> 
     
    274276         <field id="micesalt"     long_name="Mean ice salinity"                                                                                                               unit="1e-3"         /> 
    275277         <field id="miceage"      long_name="Mean ice age"                                                                                                                    unit="years"        /> 
     278         <field id="alb_ice"      long_name="Mean albedo over sea ice"                                                                                                        unit=""             /> 
     279         <field id="albedo"       long_name="Mean albedo over sea ice and ocean"                                                                                              unit=""             /> 
    276280 
    277281         <field id="iceage_cat"   long_name="Ice age for categories"                                       unit="days"   axis_ref="ncatice" /> 
     
    324328         <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   /> 
    325329         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   /> 
     330         <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(<20cm) + open water"      unit="m/day"   /> 
    326331 
    327332         <field id="afxtot"       long_name="area tendency (total)"                                        unit="day-1"   /> 
     
    365370      <field_group id="grid_U"   grid_ref="grid_U_2D"> 
    366371         <field id="e3u"          long_name="U-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_U_3D" /> 
     372         <field id="e3u_0"        long_name="Initial U-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_U_3D"/> 
    367373         <field id="utau"         long_name="Wind Stress along i-axis"                               standard_name="surface_downward_x_stress"   unit="N/m2"                            /> 
    368374         <field id="uoce"         long_name="ocean current along i-axis"                             standard_name="sea_water_x_velocity"        unit="m/s"        grid_ref="grid_U_3D" /> 
     
    400406      <field_group id="grid_V"   grid_ref="grid_V_2D"> 
    401407         <field id="e3v"          long_name="V-cell thickness"                                       standard_name="cell_thickness"              unit="m"          grid_ref="grid_V_3D" /> 
     408         <field id="e3v_0"        long_name="Initial V-cell thickness"                               standard_name="ref_cell_thickness"          unit="m"          grid_ref="grid_V_3D"/> 
    402409         <field id="vtau"         long_name="Wind Stress along j-axis"                               standard_name="surface_downward_y_stress"   unit="N/m2"                            /> 
    403410         <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" /> 
     
    454461        <!-- avt_tide: available with key_zdftmx --> 
    455462        <field id="av_tide"      long_name="tidal vertical diffusivity"   standard_name="ocean_vertical_tracer_diffusivity_due_to_tides"   unit="m2/s" /> 
     463 
     464       <!-- variables available with key_zdftmx_new --> 
     465        <field id="av_ratio"     long_name="S over T diffusivity ratio"            standard_name="salinity_over_temperature_diffusivity_ratio"                     unit="1"    /> 
     466        <field id="av_wave"      long_name="wave-induced vertical diffusivity"     standard_name="ocean_vertical_tracer_diffusivity_due_to_internal_waves"         unit="m2/s" /> 
     467        <field id="bn2"          long_name="squared Brunt-Vaisala frequency"       standard_name="squared_brunt_vaisala_frequency"                                 unit="s-1"  /> 
     468        <field id="bflx_tmx"     long_name="wave-induced buoyancy flux"            standard_name="buoyancy_flux_due_to_internal_waves"                             unit="W/kg" /> 
     469        <field id="pcmap_tmx"    long_name="power consumed by wave-driven mixing"  standard_name="vertically_integrated_power_consumption_by_wave_driven_mixing"   unit="W/m2"      grid_ref="grid_W_2D" /> 
     470        <field id="emix_tmx"     long_name="power density available for mixing"    standard_name="power_available_for_mixing_from_breaking_internal_waves"         unit="W/kg" /> 
    456471 
    457472        <!-- variables available with key_diaar5 -->    
     
    849864       <field id="Totlig"      long_name="Total ligand concentation"               unit="nmol/m3"    grid_ref="grid_T_3D" /> 
    850865       <field id="Biron"       long_name="Bioavailable iron"                       unit="nmol/m3"    grid_ref="grid_T_3D" /> 
    851        <field id="Sdenit"      long_name="Nitrate reduction in the sediments"      unit="mol/m2/s"                        /> 
     866       <field id="Sdenit"      long_name="Nitrate reduction in the sediments"      unit="molN/m2/s"                       /> 
     867       <field id="SedCal"      long_name="Calcite burial in the sediments"         unit="molC/m2/s"                       /> 
     868       <field id="SedSi"       long_name="Silicon burial in the sediments"         unit="molSi/m2/s"                      /> 
     869       <field id="SedC"        long_name="Organic C burial in the sediments"       unit="molC/m2/s"                       /> 
    852870       <field id="Ironice"     long_name="Iron input/uptake due to sea ice"        unit="mol/m2/s"                        /> 
    853871       <field id="HYDR"        long_name="Iron input from hydrothemal vents"       unit="mol/m2/s"   grid_ref="grid_T_3D" /> 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r5429 r6333  
    2121   cn_icerst_outdir = "."          !  directory in which to write output ice restarts 
    2222   ln_limdyn     = .true.          !  ice dynamics (T) or thermodynamics only (F) 
    23    rn_amax       = 0.999           !  maximum tolerated ice concentration  
     23   rn_amax_n     = 0.999           !  maximum tolerated ice concentration NH 
     24   rn_amax_s     = 0.999           !  maximum tolerated ice concentration SH 
    2425   ln_limdiahsb  = .false.         !  check the heat and salt budgets (T) or not (F) 
    2526   ln_limdiaout  = .true.          !  output the heat and salt budgets (T) or not (F) 
     
    8586   rn_hnewice  = 0.1               !  thickness for new ice formation in open water (m) 
    8687   ln_frazil   = .false.           !  use frazil ice collection thickness as a function of wind (T) or not (F) 
    87    rn_maxfrazb = 0.0               !  maximum fraction of frazil ice collecting at the ice base 
     88   rn_maxfrazb = 1.0               !  maximum fraction of frazil ice collecting at the ice base 
    8889   rn_vfrazb   = 0.417             !  thresold drift speed for frazil ice collecting at the ice bottom (m/s) 
    8990   rn_Cfrazb   = 5.0               !  squeezing coefficient for frazil ice collecting at the ice bottom 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5783 r6333  
    55!!                                    namsbc_cpl, namtra_qsr, namsbc_rnf, 
    66!!                                    namsbc_apr, namsbc_ssr, namsbc_alb) 
    7 !!              4 - lateral boundary (namlbc, namcla, namobc, namagrif, nambdy, nambdy_tide) 
     7!!              4 - lateral boundary (namlbc, namcla, namagrif, nambdy, nambdy_tide) 
    88!!              5 - bottom  boundary (nambfr, nambbc, nambbl) 
    99!!              6 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_dmp) 
    1010!!              7 - dynamics         (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf) 
    11 !!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_kpp, namzdf_ddm, namzdf_tmx) 
     11!!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_kpp, namzdf_ddm, namzdf_tmx, namzdf_tmx_new) 
    1212!!              9 - diagnostics      (namnc4, namtrd, namspr, namflo, namhsb, namsto) 
    1313!!             10 - miscellaneous    (namsol, nammpp, namctl) 
     
    408408   ln_qsr_2bd  = .false.   !  2 bands              light penetration 
    409409   ln_qsr_bio  = .false.   !  bio-model light penetration 
    410    nn_chldta   =      1    !  RGB : Chl data (=1) or cst value (=0) 
     410   nn_chldta   =      1    !  RGB : 2D Chl data (=1), 3D Chl data (=2) or cst value (=0) 
    411411   rn_abs      =   0.58    !  RGB & 2 bands: fraction of light (rn_si1) 
    412412   rn_si0      =   0.35    !  RGB & 2 bands: shortess depth of extinction 
     
    500500&namsbc_alb    !   albedo parameters 
    501501!----------------------------------------------------------------------- 
    502    rn_cloud    =    0.06   !  cloud correction to snow and ice albedo 
    503    rn_albice   =    0.53   !  albedo of melting ice in the arctic and antarctic 
    504    rn_alphd    =    0.80   !  coefficients for linear interpolation used to 
    505    rn_alphc    =    0.65   !  compute albedo between two extremes values 
    506    rn_alphdi   =    0.72   !  (Pyane, 1972) 
     502   nn_ice_alb  =    0   !  parameterization of ice/snow albedo 
     503                        !     0: Shine & Henderson-Sellers (JGR 1985) 
     504                        !     1: "home made" based on Brandt et al. (J. Climate 2005) 
     505                        !                         and Grenfell & Perovich (JGR 2004) 
     506   rn_albice   =  0.53  !  albedo of bare puddled ice (values from 0.49 to 0.58) 
     507                        !     0.53 (default) => if nn_ice_alb=0 
     508                        !     0.50 (default) => if nn_ice_alb=1 
    507509/ 
    508510!----------------------------------------------------------------------- 
     
    546548!!   namlbc        lateral momentum boundary condition 
    547549!!   namcla        cross land advection 
    548 !!   namobc        open boundaries parameters                           ("key_obc") 
    549550!!   namagrif      agrif nested grid ( read by child model only )       ("key_agrif") 
    550551!!   nambdy        Unstructured open boundaries                         ("key_bdy") 
     
    563564!----------------------------------------------------------------------- 
    564565   nn_cla      =    0      !  advection between 2 ocean pts separates by land 
    565 / 
    566 !----------------------------------------------------------------------- 
    567 &namobc        !   open boundaries parameters                           ("key_obc") 
    568 !----------------------------------------------------------------------- 
    569    ln_obc_clim = .false.   !  climatological obc data files (T) or not (F) 
    570    ln_vol_cst  = .true.    !  impose the total volume conservation (T) or not (F) 
    571    ln_obc_fla  = .false.   !  Flather open boundary condition 
    572    nn_obcdta   =    1      !  = 0 the obc data are equal to the initial state 
    573                            !  = 1 the obc data are read in 'obc.dta' files 
    574    cn_obcdta   = 'annual'  !  set to annual if obc datafile hold 1 year of data 
    575                            !  set to monthly if obc datafile hold 1 month of data 
    576    rn_dpein    =    1.     !  damping time scale for inflow at east  open boundary 
    577    rn_dpwin    =    1.     !     -           -         -       west    -      - 
    578    rn_dpnin    =    1.     !     -           -         -       north   -      - 
    579    rn_dpsin    =    1.     !     -           -         -       south   -      - 
    580    rn_dpeob    = 3000.     !  time relaxation (days) for the east  open boundary 
    581    rn_dpwob    =   15.     !     -           -         -     west    -      - 
    582    rn_dpnob    = 3000.     !     -           -         -     north   -      - 
    583    rn_dpsob    =   15.     !     -           -         -     south   -      - 
    584    rn_volemp   =    1.     !  = 0 the total volume change with the surface flux (E-P-R) 
    585                            !  = 1 the total volume remains constant 
    586566/ 
    587567!----------------------------------------------------------------------- 
     
    898878!!             Tracers & Dynamics vertical physics namelists 
    899879!!====================================================================== 
    900 !!    namzdf        vertical physics 
    901 !!    namzdf_ric    richardson number dependent vertical mixing         ("key_zdfric") 
    902 !!    namzdf_tke    TKE dependent vertical mixing                       ("key_zdftke") 
    903 !!    namzdf_kpp    KPP dependent vertical mixing                       ("key_zdfkpp") 
    904 !!    namzdf_ddm    double diffusive mixing parameterization            ("key_zdfddm") 
    905 !!    namzdf_tmx    tidal mixing parameterization                       ("key_zdftmx") 
     880!!    namzdf            vertical physics 
     881!!    namzdf_ric        richardson number dependent vertical mixing     ("key_zdfric") 
     882!!    namzdf_tke        TKE dependent vertical mixing                   ("key_zdftke") 
     883!!    namzdf_kpp        KPP dependent vertical mixing                   ("key_zdfkpp") 
     884!!    namzdf_ddm        double diffusive mixing parameterization        ("key_zdfddm") 
     885!!    namzdf_tmx        tidal mixing parameterization                   ("key_zdftmx") 
     886!!    namzdf_tmx_new    new tidal mixing parameterization               ("key_zdftmx_new") 
    906887!!====================================================================== 
    907888! 
     
    1010991   rn_tfe_itf  = 1.        !  ITF tidal dissipation efficiency 
    1011992/ 
    1012  
     993!----------------------------------------------------------------------- 
     994&namzdf_tmx_new    !   new tidal mixing parameterization                ("key_zdftmx_new") 
     995!----------------------------------------------------------------------- 
     996   nn_zpyc     = 1         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2) 
     997   ln_mevar    = .true.    !  variable (T) or constant (F) mixing efficiency 
     998   ln_tsdiff   = .true.    !  account for differential T/S mixing (T) or not (F) 
     999/ 
    10131000!!====================================================================== 
    10141001!!                  ***  Miscellaneous namelists  *** 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/CONFIG/SHARED/namelist_top_ref

    r5416 r6333  
    6262   rn_ahtrc_0       =  2000.    !  horizontal eddy diffusivity for tracers [m2/s] 
    6363   rn_ahtrb_0       =     0.    !     background eddy diffusivity for ldf_iso [m2/s] 
     64   rn_fact_lap      =     1.    !     enhanced zonal eddy diffusivity 
    6465/ 
    6566!----------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modarrays.F90

    r6237 r6333  
    11! 
    2 ! $Id: modarrays.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modbc.F90

    r6237 r6333  
    11! 
    2 ! $Id: modbc.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modbcfunction.F90

    r6237 r6333  
    11! 
    2 ! $Id: modbcfunction.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modcluster.F90

    r6237 r6333  
    11! 
    2 ! $Id: modcluster.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modcurgridfunctions.F90

    r6237 r6333  
    11! 
    2 ! $Id: modcurgridfunctions.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modinitvars.F90

    r6237 r6333  
    11! 
    2 ! $Id: modinitvars.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     Agrif (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modinterp.F90

    r6237 r6333  
    11! 
    2 ! $Id: modinterp.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modinterpbasic.F90

    r6237 r6333  
    11! 
    2 ! $Id: modinterpbasic.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modmask.F90

    r6237 r6333  
    11! 
    2 ! $Id: modmask.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modmpp.F90

    r6237 r6333  
    11! 
    2 ! $Id: modmpp.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modsauv.F90

    r6237 r6333  
    11! 
    2 ! $Id: modsauv.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modupdate.F90

    r6237 r6333  
    11! 
    2 ! $Id: modupdate.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modupdatebasic.F90

    r6237 r6333  
    11! 
    2 ! $Id: modupdatebasic.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     AGRIF (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/EXTERNAL/AGRIF/AGRIF_FILES/modutil.F90

    r6237 r6333  
    11! 
    2 ! $Id: modutil.F90 4779 2014-09-19 14:21:37Z rblod $ 
     2! $Id$ 
    33! 
    44!     Agrif (Adaptive Grid Refinement In Fortran) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90

    r5781 r6333  
    6969      IF( .NOT. ln_limini ) THEN   
    7070          
    71          tfu(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
     71         CALL eos_fzp( tsn(:,:,1,jp_sal), tfu(:,:) )       ! freezing/melting point of sea water [Celcius] 
     72         tfu(:,:) = tfu(:,:) *  tmask(:,:,1) 
    7273 
    7374         DO jj = 1, jpj 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r5781 r6333  
    301301 
    302302   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice 
     303   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_amax_2d  !: maximum ice concentration 2d array 
    303304 
    304305   !!-------------------------------------------------------------------------- 
     
    372373   INTEGER          , PUBLIC ::   nlay_i          !: number of ice  layers  
    373374   INTEGER          , PUBLIC ::   nlay_s          !: number of snow layers  
    374    CHARACTER(len=32), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input) 
     375   CHARACTER(len=80), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input) 
    375376   CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory 
    376    CHARACTER(len=32), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output) 
     377   CHARACTER(len=80), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output) 
    377378   CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory 
    378379   LOGICAL          , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F) 
    379380   LOGICAL          , PUBLIC ::   ln_icectl       !: flag for sea-ice points output (T) or not (F) 
    380    REAL(wp)         , PUBLIC ::   rn_amax         !: maximum ice concentration 
     381   REAL(wp)         , PUBLIC ::   rn_amax_n       !: maximum ice concentration Northern hemisphere 
     382   REAL(wp)         , PUBLIC ::   rn_amax_s       !: maximum ice concentration Southern hemisphere 
    381383   INTEGER          , PUBLIC ::   iiceprt         !: debug i-point 
    382384   INTEGER          , PUBLIC ::   jiceprt         !: debug j-point 
     
    438440         &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) ,                        & 
    439441         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), qlead  (jpi,jpj) ,                     & 
     442         &      rn_amax_2d(jpi,jpj),                                                            & 
    440443         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,                        & 
    441444         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,    & 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r5781 r6333  
    256256            ENDIF 
    257257            IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
    258             IF (     zamax   > rn_amax+epsi10 .AND. cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN 
     258            IF (     zamax   > MAX( rn_amax_n, rn_amax_s ) + epsi10 .AND. & 
     259               &                         cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN 
    259260                                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    260261            ENDIF 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r5781 r6333  
    117117 
    118118      ! basal temperature (considered at freezing point) 
    119       t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1)  
     119      CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
     120      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
    120121 
    121122      IF( ln_iceini ) THEN 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r5781 r6333  
    9494      !!              - fr_i    : ice fraction 
    9595      !!              - tn_ice  : sea-ice surface temperature 
    96       !!              - alb_ice : sea-ice albedo (only useful in coupled mode) 
     96      !!              - alb_ice : sea-ice albedo (recomputed only for coupled mode) 
    9797      !! 
    9898      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 
     
    106106      REAL(wp) ::   zqsr                                           ! New solar flux received by the ocean 
    107107      ! 
    108       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 2D/3D workspace 
     108      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 3D workspace 
     109      REAL(wp), POINTER, DIMENSION(:,:)   ::   zalb                 ! 2D workspace 
    109110      !!--------------------------------------------------------------------- 
    110111 
    111112      ! make calls for heat fluxes before it is modified 
     113      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
    112114      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface 
    113115      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface 
     
    118120      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
    119121         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 
    120       IF( iom_use('qemp_oce' ) )   CALL iom_put( "qemp_oce"  , qemp_oce(:,:) )   
    121       IF( iom_use('qemp_ice' ) )   CALL iom_put( "qemp_ice"  , qemp_ice(:,:) )   
    122  
    123       ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
     122      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )   
     123      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )   
     124      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   ! emp over ocean (taking into account the snow blown away from the ice) 
     125      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice) 
     126 
     127      ! clem 2016: albedo output 
     128      CALL wrk_alloc( jpi,jpj, zalb )     
     129 
     130      zalb(:,:) = 0._wp 
     131      WHERE     ( SUM( a_i_b, dim=3 ) <= epsi06 )  ;  zalb(:,:) = 0.066_wp 
     132      ELSEWHERE                                    ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
     133      END WHERE 
     134      IF( iom_use('alb_ice' ) )  CALL iom_put( "alb_ice"  , zalb(:,:) )           ! ice albedo output 
     135 
     136      zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - SUM( a_i_b, dim=3 ) )       
     137      IF( iom_use('albedo'  ) )  CALL iom_put( "albedo"  , zalb(:,:) )           ! ice albedo output 
     138 
     139      CALL wrk_dealloc( jpi,jpj, zalb )     
     140      ! 
     141       
    124142      DO jj = 1, jpj 
    125143         DO ji = 1, jpi 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r5781 r6333  
    7575      INTEGER ::   ii, ij, iter     !   -       - 
    7676      REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zde                         ! local scalars 
    77       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      - 
     77      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      - 
    7878      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
    79       LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness 
    8079      CHARACTER (len = 15) :: fieldid 
    8180 
     
    108107      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i 
    109108 
    110       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d   !: 1-D version of e_i 
    111  
    112       REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity 
    113  
    114       REAL(wp) :: zcai = 1.4e-3_wp 
     109      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d !: 1-D version of e_i 
     110 
     111      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel     ! relative ice / frazil velocity 
     112 
     113      REAL(wp) :: zcai = 1.4e-3_wp                     ! ice-air drag (clem: should be dependent on coupling/forcing used) 
    115114      !!-----------------------------------------------------------------------! 
    116115 
     
    143142      !------------------------------------------------------------------------------!     
    144143      ! hicol is the thickness of new ice formed in open water 
    145       ! hicol can be either prescribed (frazswi = 0) 
    146       ! or computed (frazswi = 1) 
     144      ! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1) 
    147145      ! Frazil ice forms in open water, is transported by wind 
    148146      ! accumulates at the edge of the consolidated ice edge 
     
    155153      zvrel(:,:) = 0._wp 
    156154 
    157       ! Default new ice thickness  
    158       hicol(:,:) = rn_hnewice 
     155      ! Default new ice thickness 
     156      WHERE( qlead(:,:) < 0._wp ) ; hicol = rn_hnewice 
     157      ELSEWHERE                   ; hicol = 0._wp 
     158      END WHERE 
    159159 
    160160      IF( ln_frazil ) THEN 
     
    182182                     &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp 
    183183                  ! Square root of wind stress 
    184                   ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) ) 
     184                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
    185185 
    186186                  !--------------------- 
     
    205205                  zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
    206206                     &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) 
    207                   zvrel(ji,jj)  = SQRT( zvrel2 ) 
     207                  zvrel(ji,jj) = SQRT( zvrel2 ) 
    208208 
    209209                  !--------------------- 
    210210                  ! Iterative procedure 
    211211                  !--------------------- 
    212                   hicol(ji,jj) = zhicrit + 0.1  
    213                   hicol(ji,jj) = zhicrit +   hicol(ji,jj)    & 
    214                      &                   / ( hicol(ji,jj) * hicol(ji,jj) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
    215  
    216 !!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1) 
    217 !!gm                                                   = zhicrit**2 + 0.2*zhicrit +0.01 
    218 !!gm                therefore the 2 lines with hicol can be replaced by 1 line: 
    219 !!gm              hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2 
    220 !!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop 
     212                  hicol(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    & 
     213                     &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
    221214 
    222215                  iter = 1 
    223                   iterate_frazil = .true. 
    224  
    225                   DO WHILE ( iter < 100 .AND. iterate_frazil )  
    226                      zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 
    227                         - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 
    228                      zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) & 
    229                         - zhicrit * ztwogp * zvrel2 
    230                      zhicol_new = hicol(ji,jj) - zf/zfp 
    231                      hicol(ji,jj)   = zhicol_new 
    232  
     216                  DO WHILE ( iter < 20 )  
     217                     zf  = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) -   & 
     218                        &    hicol(ji,jj) * zhicrit * ztwogp * zvrel2 
     219                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 
     220 
     221                     hicol(ji,jj) = hicol(ji,jj) - zf/zfp 
    233222                     iter = iter + 1 
    234  
    235                   END DO ! do while 
     223                  END DO 
    236224 
    237225               ENDIF ! end of selection of pixels where ice forms 
    238226 
    239             END DO ! loop on ji ends 
    240          END DO ! loop on jj ends 
    241       !  
    242       CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 
    243       CALL lbc_lnk( hicol(:,:), 'T', 1. ) 
     227            END DO  
     228         END DO  
     229         !  
     230         CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 
     231         CALL lbc_lnk( hicol(:,:), 'T', 1. ) 
    244232 
    245233      ENDIF ! End of computation of frazil ice collection thickness 
     
    282270      ! Move from 2-D to 1-D vectors 
    283271      !------------------------------ 
    284       ! If ocean gains heat do nothing  
    285       ! 0therwise compute new ice formation 
     272      ! If ocean gains heat do nothing. Otherwise compute new ice formation 
    286273 
    287274      IF ( nbpac > 0 ) THEN 
     
    297284         END DO 
    298285 
    299          CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) ) 
    300          CALL tab_2d_1d( nbpac, t_bo_1d   (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) ) 
    301          CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac)     , sfx_opw, jpi, jpj, npac(1:nbpac) ) 
    302          CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) ) 
    303          CALL tab_2d_1d( nbpac, hicol_1d  (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) ) 
    304          CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) ) 
    305  
    306          CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd, jpi, jpj, npac(1:nbpac) ) 
    307          CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw, jpi, jpj, npac(1:nbpac) ) 
     286         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead     , jpi, jpj, npac(1:nbpac) ) 
     287         CALL tab_2d_1d( nbpac, t_bo_1d   (1:nbpac)     , t_bo      , jpi, jpj, npac(1:nbpac) ) 
     288         CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac)     , sfx_opw   , jpi, jpj, npac(1:nbpac) ) 
     289         CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw   , jpi, jpj, npac(1:nbpac) ) 
     290         CALL tab_2d_1d( nbpac, hicol_1d  (1:nbpac)     , hicol     , jpi, jpj, npac(1:nbpac) ) 
     291         CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel     , jpi, jpj, npac(1:nbpac) ) 
     292 
     293         CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd   , jpi, jpj, npac(1:nbpac) ) 
     294         CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw   , jpi, jpj, npac(1:nbpac) ) 
     295         CALL tab_2d_1d( nbpac, rn_amax_1d(1:nbpac)     , rn_amax_2d, jpi, jpj, npac(1:nbpac) ) 
    308296 
    309297         !------------------------------------------------------------------------------! 
     
    316304         zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:)  
    317305         za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 
     306 
    318307         !---------------------- 
    319308         ! Thickness of new ice 
    320309         !---------------------- 
    321          DO ji = 1, nbpac 
    322             zh_newice(ji) = rn_hnewice 
    323          END DO 
    324          IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 
     310         zh_newice(1:nbpac) = hicol_1d(1:nbpac) 
    325311 
    326312         !---------------------- 
     
    384370            ! salt flux 
    385371            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 
    386  
     372         END DO 
     373          
     374         zv_frazb(:) = 0._wp 
     375         IF( ln_frazil ) THEN 
    387376            ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
    388             rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 
    389             zfrazb        = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 
    390             zv_frazb(ji)  =         zfrazb   * zv_newice(ji) 
    391             zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
    392          END DO 
    393  
     377            DO ji = 1, nbpac 
     378               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 
     379               zfrazb        = rswitch * ( TANH( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 
     380               zv_frazb(ji)  =         zfrazb   * zv_newice(ji) 
     381               zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
     382            END DO 
     383         END IF 
     384          
    394385         !----------------- 
    395386         ! Area of new ice 
     
    409400         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    410401         DO ji = 1, nbpac 
    411             IF ( za_newice(ji) >  ( rn_amax - zat_i_1d(ji) ) ) THEN 
    412                zda_res(ji)   = za_newice(ji) - ( rn_amax - zat_i_1d(ji) ) 
     402            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN 
     403               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) ) 
    413404               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
    414405               za_newice(ji) = za_newice(ji) - zda_res  (ji) 
     
    443434               jl = jcat(ji) 
    444435               rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 
    445                ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                      & 
     436               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    & 
    446437                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  & 
    447438                  &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 ) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5781 r6333  
    422422            DO jj = 1, jpj 
    423423               DO ji = 1, jpi 
    424                   a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax ) 
     424                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax_2d(ji,jj) ) 
    425425               END DO 
    426426            END DO 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r5781 r6333  
    8080         DO jj = 1, jpj 
    8181            DO ji = 1, jpi 
    82                IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    83                   a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    84                   oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     82               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     83                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) ) 
     84                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) ) 
    8585               ENDIF 
    8686            END DO 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r5781 r6333  
    9494         DO jj = 1, jpj 
    9595            DO ji = 1, jpi 
    96                IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    97                   a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    98                   oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     96               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     97                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) ) 
     98                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) ) 
    9999               ENDIF 
    100100            END DO 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r5781 r6333  
    157157      ENDIF 
    158158 
    159       IF ( iom_use( "icecolf" ) ) THEN  
    160          DO jj = 1, jpj 
    161             DO ji = 1, jpi 
    162                rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) ) ) 
    163                z2d(ji,jj) = hicol(ji,jj) * rswitch 
    164             END DO 
    165          END DO 
    166          CALL iom_put( "icecolf"     , z2d              )        ! frazil ice collection thickness 
    167       ENDIF 
     159      IF ( iom_use( "icecolf" ) )   CALL iom_put( "icecolf", hicol )  ! frazil ice collection thickness 
    168160 
    169161      CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature 
     
    235227      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice  
    236228      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
     229 
     230 
     231      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
     232         DO jj = 1, jpj  
     233            DO ji = 1, jpi 
     234               z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) ! mean ice thickness 
     235            END DO 
     236         END DO 
     237         WHERE( z2d(:,:) < 0.2 .AND. z2d(:,:) > 0. ) ; z2da = wfx_bog 
     238         ELSEWHERE                                   ; z2da = 0._wp 
     239         END WHERE 
     240         CALL iom_put( "vfxthin", ( wfx_opw + z2da ) * ztmp ) 
     241      ENDIF 
    237242       
    238243      !-------------------------------- 
     
    311316      !! 
    312317      !! History : 
    313       !!   4.1  !  2013-06  (C. Rousset) 
     318      !!   4.0  !  2013-06  (C. Rousset) 
    314319      !!---------------------------------------------------------------------- 
    315320      INTEGER, INTENT( in ) ::   kt               ! ocean time-step index) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r5781 r6333  
    5151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_dif_1d 
    5252   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_opw_1d 
     53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d 
    5354   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d 
    5455   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_1d 
     
    144145         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,    &  
    145146         &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) ,                      & 
     147         &      rn_amax_1d(jpij) ,                                         & 
    146148         &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      & 
    147149         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) ,   & 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r5781 r6333  
    658658 
    659659      DO jk = 1, jpkm1 
    660          fzptnz(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) ) 
     660        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) ) 
    661661      END DO 
    662662 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r6237 r6333  
    145145      ENDIF 
    146146 
     147      ! Output of initial vertical scale factor 
     148      CALL iom_put("e3t_0", e3t_0(:,:,:) ) 
     149      CALL iom_put("e3u_0", e3t_0(:,:,:) ) 
     150      CALL iom_put("e3v_0", e3t_0(:,:,:) ) 
     151      ! 
    147152      IF( .NOT.lk_vvl ) THEN 
    148153         CALL iom_put( "e3t" , fse3t_n(:,:,:) ) 
     
    153158 
    154159      CALL iom_put( "ssh" , sshn )                 ! sea surface height 
    155       if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
    156160       
    157161      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r6237 r6333  
    139139      ! horizontal grid definition 
    140140 
    141 #if ! defined key_xios2 
    142141      CALL set_scalar 
    143 #endif 
    144142 
    145143      IF( TRIM(cdname) == TRIM(cxios_context) ) THEN   
     
    11931191      REAL(wp), DIMENSION(:)   , OPTIONAL, INTENT(in) ::   lonvalue, latvalue 
    11941192      REAL(wp), DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   bounds_lon, bounds_lat, area 
    1195       LOGICAL,  DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   mask 
     1193#if ! defined key_xios2 
     1194     LOGICAL,  DIMENSION(:,:) , OPTIONAL, INTENT(in) ::   mask 
     1195#else 
     1196      LOGICAL,  DIMENSION(:) , OPTIONAL, INTENT(in) ::   mask 
     1197#endif 
    11961198 
    11971199#if ! defined key_xios2 
     
    12151217         CALL xios_set_domain_attr     ( cdid, ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj,   & 
    12161218            &    data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj ,   & 
    1217             &    lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_2D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon,      & 
     1219            &    lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_1D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon,                  & 
    12181220            &    bounds_lat_1D=bounds_lat, area=area, type='curvilinear') 
    12191221     ENDIF 
     
    12211223         CALL xios_set_domaingroup_attr( cdid, ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj,   & 
    12221224            &    data_dim=data_dim, data_ibegin=data_ibegin, data_ni=data_ni, data_jbegin=data_jbegin, data_nj=data_nj ,   & 
    1223             &    lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_2D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon,      & 
     1225            &    lonvalue_1D=lonvalue, latvalue_1D=latvalue, mask_1D=mask, nvertex=nvertex, bounds_lon_1D=bounds_lon,                  & 
    12241226            &    bounds_lat_1D=bounds_lat, area=area, type='curvilinear' ) 
    12251227      ENDIF 
     
    12341236     INTEGER                  , OPTIONAL, INTENT(in) ::   ibegin, jbegin, ni, nj 
    12351237 
    1236      IF ( xios_is_valid_domain     (cdid) ) THEN 
     1238     IF ( xios_is_valid_zoom_domain     (cdid) ) THEN 
    12371239         CALL xios_set_zoom_domain_attr     ( cdid, ibegin=ibegin, jbegin=jbegin, ni=ni,    & 
    12381240           &   nj=nj) 
     
    13261328      IF ( xios_is_valid_gridgroup(cdid) )   CALL xios_set_gridgroup_attr( cdid, mask=mask ) 
    13271329#else 
    1328       IF ( xios_is_valid_grid     (cdid) )   CALL xios_set_grid_attr     ( cdid, mask3=mask ) 
    1329       IF ( xios_is_valid_gridgroup(cdid) )   CALL xios_set_gridgroup_attr( cdid, mask3=mask ) 
     1330      IF ( xios_is_valid_grid     (cdid) )   CALL xios_set_grid_attr     ( cdid, mask_3D=mask ) 
     1331      IF ( xios_is_valid_gridgroup(cdid) )   CALL xios_set_gridgroup_attr( cdid, mask_3D=mask ) 
    13301332#endif 
    13311333      CALL xios_solve_inheritance() 
     
    13881390         END SELECT 
    13891391         ! 
     1392#if ! defined key_xios2 
    13901393         CALL iom_set_domain_attr( "grid_"//cdgrd       , mask = RESHAPE(zmask(nldi:nlei,nldj:nlej,1),(/ni,nj    /)) /= 0. ) 
     1394#else 
     1395         CALL iom_set_domain_attr( "grid_"//cdgrd       , mask = RESHAPE(zmask(nldi:nlei,nldj:nlej,1),(/ni*nj    /)) /= 0. ) 
     1396#endif   
    13911397         CALL iom_set_grid_attr  ( "grid_"//cdgrd//"_3D", mask = RESHAPE(zmask(nldi:nlei,nldj:nlej,:),(/ni,nj,jpk/)) /= 0. ) 
    13921398      ENDIF 
     
    15321538#else 
    15331539! Pas teste : attention aux indices ! 
    1534       CALL iom_set_domain_attr("ptr", ni_glo=jpiglo, nj_glo=jpjglo, ibegin=nimpp+nldi-2, jbegin=njmpp+nldj-2, ni=ni, nj=nj) 
    1535       CALL iom_set_domain_attr("ptr", data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj) 
    1536       CALL iom_set_domain_attr("ptr", lonvalue = zlon,   & 
     1540      CALL iom_set_domain_attr("gznl", ni_glo=jpiglo, nj_glo=jpjglo, ibegin=nimpp+nldi-2, jbegin=njmpp+nldj-2, ni=ni, nj=nj) 
     1541      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = 1-nldi, data_ni = jpi, data_jbegin = 1-nldj, data_nj = jpj) 
     1542      CALL iom_set_domain_attr("gznl", lonvalue = zlon,   & 
    15371543         &                             latvalue = RESHAPE(plat(nldi:nlei, nldj:nlej),(/ ni*nj /)))   
    1538        CALL iom_set_zoom_domain_attr ('ptr', ibegin=ix, nj=jpjglo) 
     1544       CALL iom_set_zoom_domain_attr ("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=jpjglo) 
    15391545#endif 
    15401546 
     
    15521558      REAL(wp), DIMENSION(1)   ::   zz = 1. 
    15531559      !!---------------------------------------------------------------------- 
     1560#if ! defined key_xios2 
    15541561      CALL iom_set_domain_attr('scalarpoint', ni_glo=jpnij, nj_glo=1, ibegin=narea, jbegin=1, ni=1, nj=1) 
     1562#else 
     1563      CALL iom_set_domain_attr('scalarpoint', ni_glo=jpnij, nj_glo=1, ibegin=narea-1, jbegin=0, ni=1, nj=1) 
     1564#endif 
    15551565      CALL iom_set_domain_attr('scalarpoint', data_dim=2, data_ibegin = 1, data_ni = 1, data_jbegin = 1, data_nj = 1) 
    15561566       
     
    17781788            idx = INDEX(clname,'@freq@') + INDEX(clname,'@FREQ@') 
    17791789            DO WHILE ( idx /= 0 )  
    1780               IF ( output_freq%hour /= 0 ) THEN 
     1790              IF ( output_freq%timestep /= 0) THEN 
     1791                  WRITE(clfreq,'(I18,A2)')INT(output_freq%timestep),'ts'  
     1792                  itrlen = LEN_TRIM(ADJUSTL(clfreq)) 
     1793              ELSE IF ( output_freq%hour /= 0 ) THEN 
    17811794                  WRITE(clfreq,'(I19,A1)')INT(output_freq%hour),'h'  
    17821795                  itrlen = LEN_TRIM(ADJUSTL(clfreq)) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90

    r5783 r6333  
    4141 
    4242   REAL(wp), PUBLIC ::   rldf                        !: multiplicative factor of diffusive coefficient 
     43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   r_fact_lap 
    4344                                                     !: Needed to define the ratio between passive and active tracer diffusion coef.  
    4445 
     
    9293      !!                 ***  FUNCTION ldftra_oce_alloc  *** 
    9394     !!---------------------------------------------------------------------- 
    94      INTEGER, DIMENSION(3) :: ierr 
     95     INTEGER, DIMENSION(4) :: ierr 
    9596     !!---------------------------------------------------------------------- 
    9697     ierr(:) = 0 
     
    116117# endif 
    117118#endif 
     119      ALLOCATE( r_fact_lap(jpi,jpj,jpk), STAT=ierr(4) ) 
    118120      ldftra_oce_alloc = MAXVAL( ierr ) 
    119121      IF( ldftra_oce_alloc /= 0 )   CALL ctl_warn('ldftra_oce_alloc: failed to allocate arrays') 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_substitute.h90

    r5783 r6333  
    1313!   'key_traldf_c3d' :                 aht: 3D coefficient 
    1414#       define   fsahtt(i,j,k)   rldf * ahtt(i,j,k) 
    15 #       define   fsahtu(i,j,k)   rldf * ahtu(i,j,k) 
     15#       define   fsahtu(i,j,k)   rldf * ahtu(i,j,k) * r_fact_lap(i,j,k) 
    1616#       define   fsahtv(i,j,k)   rldf * ahtv(i,j,k) 
    1717#       define   fsahtw(i,j,k)   rldf * ahtw(i,j,k) 
     
    1919!   'key_traldf_c2d' :                 aht: 2D coefficient 
    2020#       define   fsahtt(i,j,k)   rldf * ahtt(i,j) 
    21 #       define   fsahtu(i,j,k)   rldf * ahtu(i,j) 
     21#       define   fsahtu(i,j,k)   rldf * ahtu(i,j) * r_fact_lap(i,j,k) 
    2222#       define   fsahtv(i,j,k)   rldf * ahtv(i,j) 
    2323#       define   fsahtw(i,j,k)   rldf * ahtw(i,j) 
     
    2525!   'key_traldf_c1d' :                aht: 1D coefficient 
    2626#       define   fsahtt(i,j,k)   rldf * ahtt(k) 
    27 #       define   fsahtu(i,j,k)   rldf * ahtu(k) 
     27#       define   fsahtu(i,j,k)   rldf * ahtu(k) * r_fact_lap(i,j,k) 
    2828#       define   fsahtv(i,j,k)   rldf * ahtv(k) 
    2929#       define   fsahtw(i,j,k)   rldf * ahtw(k) 
     
    3131!   Default option :             aht: Constant coefficient 
    3232#      define   fsahtt(i,j,k)   rldf * aht0 
    33 #      define   fsahtu(i,j,k)   rldf * aht0 
     33#      define   fsahtu(i,j,k)   rldf * aht0 * r_fact_lap(i,j,k) 
    3434#      define   fsahtv(i,j,k)   rldf * aht0 
    3535#      define   fsahtw(i,j,k)   rldf * aht0 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

    r5781 r6333  
    99   !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0 
    1010   !!             -   ! 2006-08  (G. Madec)  cleaning for surface module 
     11   !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo 
    1112   !!---------------------------------------------------------------------- 
    1213 
     
    2930 
    3031   INTEGER  ::   albd_init = 0      !: control flag for initialization 
    31    REAL(wp) ::   zzero     = 0.e0   ! constant values 
    32    REAL(wp) ::   zone      = 1.e0   !    "       " 
    33  
    34    REAL(wp) ::   c1     = 0.05    ! constants values 
    35    REAL(wp) ::   c2     = 0.10    !    "        " 
    36    REAL(wp) ::   rmue   = 0.40    !  cosine of local solar altitude 
    37  
     32   
     33   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude 
     34   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
     35   REAL(wp) ::   c1       = 0.05    ! snow thickness (only for nn_ice_alb=0) 
     36   REAL(wp) ::   c2       = 0.10    !  "        " 
     37   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0) 
     38  
    3839   !                             !!* namelist namsbc_alb 
    39    REAL(wp) ::   rn_cloud         !  cloudiness effect on snow or ice albedo (Grenfell & Perovich, 1984) 
    40 #if defined key_lim3 
    41    REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers) 
    42 #else 
    43    REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers) 
    44 #endif 
    45    REAL(wp) ::   rn_alphd         !  coefficients for linear interpolation used to compute 
    46    REAL(wp) ::   rn_alphdi        !  albedo between two extremes values (Pyane, 1972) 
    47    REAL(wp) ::   rn_alphc         !  
     40   INTEGER  ::   nn_ice_alb 
     41   REAL(wp) ::   rn_albice 
    4842 
    4943   !!---------------------------------------------------------------------- 
     
    5953      !!           
    6054      !! ** Purpose :   Computation of the albedo of the snow/ice system  
    61       !!                as well as the ocean one 
    6255      !!        
    63       !! ** Method  : - Computation of the albedo of snow or ice (choose the  
    64       !!                rignt one by a large number of tests 
    65       !!              - Computation of the albedo of the ocean 
    66       !! 
    67       !! References :   Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250. 
     56      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb) 
     57      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies 
     58      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005) 
     59      !!                                                                           and Grenfell & Perovich (JGR 2004) 
     60      !!                Description of scheme 1: 
     61      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005) 
     62      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999 
     63      !!                     0-5cm  : linear function of ice thickness 
     64      !!                     5-150cm: log    function of ice thickness 
     65      !!                     > 150cm: constant 
     66      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004) 
     67      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting 
     68      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004) 
     69      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law 
     70      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice 
     71      !! 
     72      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions: 
     73      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo 
     74      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger  
     75      !!                     under melting conditions than under freezing conditions 
     76      !!                  3) the evolution of ice albedo as a function of ice thickness shows   
     77      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic 
     78      !! 
     79      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250. 
     80      !!                Brandt et al. 2005, J. Climate, vol 18 
     81      !!                Grenfell & Perovich 2004, JGR, vol 109  
    6882      !!---------------------------------------------------------------------- 
    6983      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin) 
     
    7387      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky 
    7488      !! 
    75       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    76       INTEGER  ::   ijpl          ! number of ice categories (3rd dim of ice input arrays) 
    77       REAL(wp) ::   zalbpsnm      ! albedo of ice under clear sky when snow is melting 
    78       REAL(wp) ::   zalbpsnf      ! albedo of ice under clear sky when snow is freezing 
    79       REAL(wp) ::   zalbpsn       ! albedo of snow/ice system when ice is coverd by snow 
    80       REAL(wp) ::   zalbpic       ! albedo of snow/ice system when ice is free of snow 
    81       REAL(wp) ::   zithsn        ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow) 
    82       REAL(wp) ::   zitmlsn       ! = 1 freezinz snow (pt_ice >=rt0_snow) ; = 0 melting snow (pt_ice<rt0_snow) 
    83       REAL(wp) ::   zihsc1        ! = 1 hsn <= c1 ; = 0 hsn > c1 
    84       REAL(wp) ::   zihsc2        ! = 1 hsn >= c2 ; = 0 hsn < c2 
    85       !! 
    86       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalbfz    ! = rn_alphdi for freezing ice ; = rn_albice for melting ice 
    87       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zficeth   !  function of ice thickness 
     89      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
     90      INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays) 
     91      REAL(wp)            ::   ralb_im, ralb_sf, ralb_sm, ralb_if 
     92      REAL(wp)            ::   zswitch, z1_c1, z1_c2 
     93      REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total 
     94      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free) 
    8895      !!--------------------------------------------------------------------- 
    89        
     96 
    9097      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories 
    91  
    92       CALL wrk_alloc( jpi,jpj,ijpl, zalbfz, zficeth ) 
     98       
     99      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it ) 
    93100 
    94101      IF( albd_init == 0 )   CALL albedo_init      ! initialization  
    95102 
    96       !--------------------------- 
    97       !  Computation of  zficeth 
    98       !--------------------------- 
    99       ! ice free of snow and melts 
    100       WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice 
    101       ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi 
    102       END  WHERE 
    103  
    104       WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz 
    105       ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 ) 
    106       ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              & 
    107          &                                                                 - 0.8608 * ph_ice * ph_ice     & 
    108          &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice 
    109       ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice 
    110       END WHERE 
    111  
    112 !!gm old code 
    113 !      DO jl = 1, ijpl 
    114 !         DO jj = 1, jpj 
    115 !            DO ji = 1, jpi 
    116 !               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN 
    117 !                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl) 
    118 !               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN 
    119 !                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 ) 
    120 !               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN 
    121 !                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               & 
    122 !                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 & 
    123 !                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl) 
    124 !               ELSE 
    125 !                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl)  
    126 !               ENDIF 
    127 !            END DO 
    128 !         END DO 
    129 !      END DO 
    130 !!gm end old code 
    131        
    132       !-----------------------------------------------  
    133       !    Computation of the snow/ice albedo system  
    134       !-------------------------- --------------------- 
    135        
    136       !    Albedo of snow-ice for clear sky. 
    137       !-----------------------------------------------     
    138       DO jl = 1, ijpl 
    139          DO jj = 1, jpj 
    140             DO ji = 1, jpi 
    141                !  Case of ice covered by snow.              
    142                !                                        !  freezing snow         
    143                zihsc1   = 1.0 - MAX( zzero , SIGN( zone , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
    144                zalbpsnf = ( 1.0 - zihsc1 ) * (  zficeth(ji,jj,jl)                                             & 
    145                   &                           + ph_snw(ji,jj,jl) * ( rn_alphd - zficeth(ji,jj,jl) ) / c1  )   & 
    146                   &     +         zihsc1   * rn_alphd   
    147                !                                        !  melting snow                 
    148                zihsc2   = MAX( zzero , SIGN( zone , ph_snw(ji,jj,jl) - c2 ) ) 
    149                zalbpsnm = ( 1.0 - zihsc2 ) * ( rn_albice + ph_snw(ji,jj,jl) * ( rn_alphc - rn_albice ) / c2 )   & 
    150                   &     +         zihsc2   *   rn_alphc  
    151                ! 
    152                zitmlsn  =  MAX( zzero , SIGN( zone , pt_ice(ji,jj,jl) - rt0_snow ) )    
    153                zalbpsn  =  zitmlsn * zalbpsnm + ( 1.0 - zitmlsn ) * zalbpsnf 
    154              
    155                !  Case of ice free of snow. 
    156                zalbpic  = zficeth(ji,jj,jl)  
    157              
    158                ! albedo of the system    
    159                zithsn   = 1.0 - MAX( zzero , SIGN( zone , - ph_snw(ji,jj,jl) ) ) 
    160                pa_ice_cs(ji,jj,jl) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic 
     103       
     104      SELECT CASE ( nn_ice_alb ) 
     105 
     106      !------------------------------------------ 
     107      !  Shine and Henderson-Sellers (1985) 
     108      !------------------------------------------ 
     109      CASE( 0 ) 
     110        
     111         ralb_sf = 0.80       ! dry snow 
     112         ralb_sm = 0.65       ! melting snow 
     113         ralb_if = 0.72       ! bare frozen ice 
     114         ralb_im = rn_albice  ! bare puddled ice  
     115          
     116         !  Computation of ice albedo (free of snow) 
     117         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im 
     118         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if 
     119         END  WHERE 
     120       
     121         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
     122         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
     123         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              & 
     124            &                                                                 - 0.8608 * ph_ice * ph_ice     & 
     125            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice 
     126         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice 
     127         END WHERE 
     128      
     129         DO jl = 1, ijpl 
     130            DO jj = 1, jpj 
     131               DO ji = 1, jpi 
     132                  ! freezing snow 
     133                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 
     134                  !                                        !  freezing snow         
     135                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
     136                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
     137                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   & 
     138                     &        +         zswitch   * ralb_sf   
     139 
     140                  ! melting snow 
     141                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2 
     142                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 
     143                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   & 
     144                      &     +         zswitch   *   ralb_sm  
     145                  ! 
     146                  ! snow albedo 
     147                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     148                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
     149                
     150                  ! Ice/snow albedo 
     151                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
     152                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl) 
     153                  ! 
     154               END DO 
    161155            END DO 
    162156         END DO 
    163       END DO 
    164        
    165       !    Albedo of snow-ice for overcast sky. 
    166       !----------------------------------------------   
    167       pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rn_cloud       ! Oberhuber correction 
    168       ! 
    169       CALL wrk_dealloc( jpi,jpj,ijpl, zalbfz, zficeth ) 
     157 
     158         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky 
     159 
     160      !------------------------------------------ 
     161      !  New parameterization (2016) 
     162      !------------------------------------------ 
     163      CASE( 1 )  
     164 
     165         ralb_im = rn_albice  ! bare puddled ice 
     166! compilation of values from literature 
     167         ralb_sf = 0.85      ! dry snow 
     168         ralb_sm = 0.75      ! melting snow 
     169         ralb_if = 0.60      ! bare frozen ice 
     170! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 
     171!         ralb_sf = 0.85       ! dry snow 
     172!         ralb_sm = 0.72       ! melting snow 
     173!         ralb_if = 0.65       ! bare frozen ice 
     174! Brandt et al 2005 (East Antarctica) 
     175!         ralb_sf = 0.87      ! dry snow 
     176!         ralb_sm = 0.82      ! melting snow 
     177!         ralb_if = 0.54      ! bare frozen ice 
     178!  
     179         !  Computation of ice albedo (free of snow) 
     180         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
     181         z1_c2 = 1. / 0.05 
     182         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im 
     183         ELSE WHERE                                              ;   zalb = ralb_if 
     184         END  WHERE 
     185          
     186         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
     187         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  & 
     188            &                                                                     ( LOG(1.5) - LOG(ph_ice) ) 
     189         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice 
     190         END WHERE 
     191 
     192         z1_c1 = 1. / 0.02 
     193         z1_c2 = 1. / 0.03 
     194         !  Computation of the snow/ice albedo 
     195         DO jl = 1, ijpl 
     196            DO jj = 1, jpj 
     197               DO ji = 1, jpi 
     198                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 ); 
     199                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 ); 
     200 
     201                   ! snow albedo 
     202                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     203                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
     204 
     205                  ! Ice/snow albedo    
     206                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
     207                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl) 
     208 
     209              END DO 
     210            END DO 
     211         END DO 
     212         ! Effect of the clouds (2d order polynomial) 
     213         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 );  
     214 
     215      END SELECT 
     216       
     217      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it ) 
    170218      ! 
    171219   END SUBROUTINE albedo_ice 
     
    181229      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky 
    182230      !! 
    183       REAL(wp) ::   zcoef   ! local scalar 
    184       !!---------------------------------------------------------------------- 
    185       ! 
    186       zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )      ! Parameterization of Briegled and Ramanathan, 1982  
    187       pa_oce_cs(:,:) = zcoef                
    188       pa_oce_os(:,:)  = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972 
     231      REAL(wp) :: zcoef  
     232      !!---------------------------------------------------------------------- 
     233      ! 
     234      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982 
     235      pa_oce_cs(:,:) = zcoef  
     236      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972 
    189237      ! 
    190238   END SUBROUTINE albedo_oce 
     
    200248      !!---------------------------------------------------------------------- 
    201249      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    202       NAMELIST/namsbc_alb/ rn_cloud, rn_albice, rn_alphd, rn_alphdi, rn_alphc 
     250      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice  
    203251      !!---------------------------------------------------------------------- 
    204252      ! 
     
    219267         WRITE(numout,*) '~~~~~~~' 
    220268         WRITE(numout,*) '   Namelist namsbc_alb : albedo ' 
    221          WRITE(numout,*) '      correction for snow and ice albedo                  rn_cloud  = ', rn_cloud 
    222          WRITE(numout,*) '      albedo of melting ice in the arctic and antarctic   rn_albice = ', rn_albice 
    223          WRITE(numout,*) '      coefficients for linear                             rn_alphd  = ', rn_alphd 
    224          WRITE(numout,*) '      interpolation used to compute albedo                rn_alphdi = ', rn_alphdi 
    225          WRITE(numout,*) '      between two extremes values (Pyane, 1972)           rn_alphc  = ', rn_alphc 
     269         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb 
     270         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice 
    226271      ENDIF 
    227272      ! 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90

    r5781 r6333  
    103103                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 ) 
    104104          
    105          fr_i(:,:) = eos_fzp( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius] 
     105         CALL eos_fzp( sss_m(:,:), fr_i(:,:) )       ! sea surface freezing temperature [Celcius] 
     106         fr_i(:,:) = fr_i(:,:) * tmask(:,:,1) 
    106107 
    107108         IF( ln_cpl )   a_i(:,:,1) = fr_i(:,:)          
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r5781 r6333  
    110110      INTEGER  ::   jl                 ! dummy loop index 
    111111      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
    112       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled) 
    113112      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice  
    114113      !!---------------------------------------------------------------------- 
     
    126125          
    127126         ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
    128          t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )   
    129           
     127         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) 
     128         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 
     129           
    130130         ! Mask sea ice surface temperature (set to rt0 over land) 
    131131         DO jl = 1, jpl 
     
    196196         ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%] 
    197197         !---------------------------------------------------------------------------------------- 
    198          CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     198         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs ) 
    199199         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    200200 
     
    202202         CASE( jp_clio )                                       ! CLIO bulk formulation 
    203203            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
    204             ! (zalb_ice) is computed within the bulk routine 
    205             CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice ) 
    206             IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
    207             IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     204            ! (alb_ice) is computed within the bulk routine 
     205            CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice ) 
     206            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     207            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    208208         CASE( jp_core )                                       ! CORE bulk formulation 
    209209            ! albedo depends on cloud fraction because of non-linear spectral effects 
    210             zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    211             CALL blk_ice_core_flx( t_su, zalb_ice ) 
    212             IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
    213             IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     210            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     211            CALL blk_ice_core_flx( t_su, alb_ice ) 
     212            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     213            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    214214         CASE ( jp_purecpl ) 
    215215            ! albedo depends on cloud fraction because of non-linear spectral effects 
    216             zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    217                                  CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su ) 
     216            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     217                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su ) 
    218218            ! clem: evap_ice is forced to 0 in coupled mode for now  
    219219            !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models 
    220220            evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp 
    221             IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
     221            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx ) 
    222222         END SELECT 
    223          CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     223         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs ) 
    224224 
    225225         !----------------------------! 
     
    264264      !!---------------------------------------------------------------------- 
    265265      INTEGER :: ierr 
     266      INTEGER :: ji, jj 
    266267      !!---------------------------------------------------------------------- 
    267268      IF(lwp) WRITE(numout,*) 
     
    320321      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
    321322      ! 
     323      DO jj = 1, jpj 
     324         DO ji = 1, jpi 
     325            IF( gphit(ji,jj) > 0._wp ) THEN  ;  rn_amax_2d(ji,jj) = rn_amax_n  ! NH 
     326            ELSE                             ;  rn_amax_2d(ji,jj) = rn_amax_s  ! SH 
     327            ENDIF 
     328        ENDDO 
     329      ENDDO  
     330      ! 
    322331      nstart = numit  + nn_fsbc       
    323332      nitrun = nitend - nit000 + 1  
     
    342351      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    343352      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  & 
    344          &                ln_limdyn, rn_amax, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt   
     353         &                ln_limdyn, rn_amax_n, rn_amax_s, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt   
    345354      !!------------------------------------------------------------------- 
    346355      !                     
     
    363372         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s 
    364373         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
    365          WRITE(numout,*) '   maximum ice concentration                               = ', rn_amax  
     374         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n  
     375         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s 
    366376         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb 
    367377         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90

    r5781 r6333  
    150150 
    151151         ! ... masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
    152          tfu(:,:) = eos_fzp( sss_m ) +  rt0  
     152         CALL eos_fzp( sss_m(:,:), tfu(:,:) ) 
     153         tfu(:,:) = tfu(:,:) + rt0 
    153154 
    154155         zsist (:,:,1) = sist (:,:) + rt0 * ( 1. - tmask(:,:,1) ) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r6237 r6333  
    370370             ! Calculate freezing temperature 
    371371                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04  
    372                 zt_frz = eos_fzp(tsb(ji,jj,ik,jp_sal), zpress)  
     372                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress)  
    373373                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp 
    374374             ENDDO 
     
    452452      zti(:,:)=tinsitu( ttbl, stbl, zpress ) 
    453453! Calculate freezing temperature 
    454       zfrz(:,:)=eos_fzp( sss_m(:,:), zpress ) 
     454      CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress ) 
    455455 
    456456       
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r6237 r6333  
    5252   REAL(wp)                   ::   rn_hrnf         !: runoffs, depth over which enhanced vertical mixing is used 
    5353   REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s] 
    54    REAL(wp)                  ::   rn_rfact        !: multiplicative factor for runoff 
     54   REAL(wp)          , PUBLIC ::   rn_rfact        !: multiplicative factor for runoff 
    5555 
    5656   LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/SOL/solver.F90

    r5783 r6333  
    9292      IF( .NOT. lk_agrif .OR. .NOT. ln_rstart) THEN 
    9393         IF( sol_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'solver_init : unable to allocate sol_oce arrays' ) 
     94         gcx (:,:) = 0.e0 
     95         gcxb(:,:) = 0.e0 
    9496      ENDIF 
    9597 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r5781 r6333  
    2222   !!             -   ! 2013-04  (F. Roquet, G. Madec)  add eos_rab, change bn2 computation and reorganize the module 
    2323   !!             -   ! 2014-09  (F. Roquet)  add TEOS-10, S-EOS, and modify EOS-80 
     24   !!             -   ! 2015-06  (P.A. Bouttier) eos_fzp functions changed to subroutines for AGRIF 
    2425   !!---------------------------------------------------------------------- 
    2526 
     
    991992 
    992993 
    993    FUNCTION eos_fzp_2d( psal, pdep ) RESULT( ptf ) 
     994   SUBROUTINE  eos_fzp_2d( psal, ptf, pdep ) 
    994995      !!---------------------------------------------------------------------- 
    995996      !!                 ***  ROUTINE eos_fzp  *** 
     
    10051006      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu] 
    10061007      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m] 
    1007       REAL(wp), DIMENSION(jpi,jpj)                          ::   ptf   ! freezing temperature [Celcius] 
     1008      REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celcius] 
    10081009      ! 
    10091010      INTEGER  ::   ji, jj   ! dummy loop indices 
     
    10381039         nstop = nstop + 1 
    10391040         ! 
    1040       END SELECT 
    1041       ! 
    1042    END FUNCTION eos_fzp_2d 
    1043  
    1044   FUNCTION eos_fzp_0d( psal, pdep ) RESULT( ptf ) 
     1041      END SELECT       
     1042      ! 
     1043  END SUBROUTINE eos_fzp_2d 
     1044 
     1045  SUBROUTINE eos_fzp_0d( psal, ptf, pdep ) 
    10451046      !!---------------------------------------------------------------------- 
    10461047      !!                 ***  ROUTINE eos_fzp  *** 
     
    10541055      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978 
    10551056      !!---------------------------------------------------------------------- 
    1056       REAL(wp), INTENT(in)           ::   psal   ! salinity   [psu] 
    1057       REAL(wp), INTENT(in), OPTIONAL ::   pdep   ! depth      [m] 
    1058       REAL(wp)                       ::   ptf   ! freezing temperature [Celcius] 
     1057      REAL(wp), INTENT(in )           ::   psal         ! salinity   [psu] 
     1058      REAL(wp), INTENT(in ), OPTIONAL ::   pdep         ! depth      [m] 
     1059      REAL(wp), INTENT(out)           ::   ptf          ! freezing temperature [Celcius] 
    10591060      ! 
    10601061      REAL(wp) :: zs   ! local scalars 
     
    10861087      END SELECT 
    10871088      ! 
    1088    END FUNCTION eos_fzp_0d 
     1089   END SUBROUTINE eos_fzp_0d 
    10891090 
    10901091 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90

    r5781 r6333  
    173173         END DO  
    174174      END DO  
    175       zfzp(:,:) = eos_fzp( tsn(:,:,1,jp_sal), zpres(:,:) ) 
     175      CALL eos_fzp( tsn(:,:,1,jp_sal), zfzp(:,:), zpres(:,:) ) 
    176176      DO jk = 1, jpk 
    177177         DO jj = 1, jpj 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r5783 r6333  
    6868      ! 
    6969      rldf = 1     ! For active tracers the  
     70      r_fact_lap(:,:,:) = 1.0 
    7071 
    7172      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r5781 r6333  
    1010   !!             -   !  2005-11  (G. Madec) zco, zps, sco coordinate 
    1111   !!            3.2  !  2009-04  (G. Madec & NEMO team)  
    12    !!            4.0  !  2012-05  (C. Rousset) store attenuation coef for use in ice model  
     12   !!            3.4  !  2012-05  (C. Rousset) store attenuation coef for use in ice model  
     13   !!            3.6  !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 
    1314   !!---------------------------------------------------------------------- 
    1415 
     
    9394      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 
    9495      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 
     96      !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562 
    9597      !!---------------------------------------------------------------------- 
    9698      ! 
     
    101103      REAL(wp) ::   zchl, zcoef, zfact   ! local scalars 
    102104      REAL(wp) ::   zc0, zc1, zc2, zc3   !    -         - 
    103       REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         - 
    104105      REAL(wp) ::   zz0, zz1, z1_e3t     !    -         - 
     106      REAL(wp) ::   zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze 
     107      REAL(wp) ::   zlogc, zlogc2, zlogc3  
    105108      REAL(wp), POINTER, DIMENSION(:,:  ) :: zekb, zekg, zekr 
    106       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 
    107       !!---------------------------------------------------------------------- 
     109      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt, zchl3d 
     110      !!-------------------------------------------------------------------------- 
    108111      ! 
    109112      IF( nn_timing == 1 )  CALL timing_start('tra_qsr') 
    110113      ! 
    111114      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    112       CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
     115      CALL wrk_alloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )  
    113116      ! 
    114117      IF( kt == nit000 ) THEN 
     
    183186            !                                             ! ------------------------- ! 
    184187            ! Set chlorophyl concentration 
    185             IF( nn_chldta == 1 .OR. lk_vvl ) THEN            !*  Variable Chlorophyll or ocean volume 
    186                ! 
    187                IF( nn_chldta == 1 ) THEN                             !*  Variable Chlorophyll 
    188                   ! 
    189                   CALL fld_read( kt, 1, sf_chl )                         ! Read Chl data and provides it at the current time step 
    190                   !          
    191 !CDIR COLLAPSE 
     188            IF( nn_chldta == 1 .OR. nn_chldta == 2 .OR. lk_vvl ) THEN    !*  Variable Chlorophyll or ocean volume 
     189               ! 
     190               IF( nn_chldta == 1 ) THEN        !*  2D Variable Chlorophyll 
     191                  ! 
     192                  CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step 
     193                  DO jk = 1, nksr + 1 
     194                     zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1)  
     195                  ENDDO 
     196                  ! 
     197               ELSE IF( nn_chldta == 2 ) THEN    !*   -3-D Variable Chlorophyll 
     198                  ! 
     199                  CALL fld_read( kt, 1, sf_chl )            ! Read Chl data and provides it at the current time step 
     200!CDIR NOVERRCHK   ! 
     201                  DO jj = 1, jpj 
    192202!CDIR NOVERRCHK 
    193                   DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl 
    194 !CDIR NOVERRCHK 
    195                      DO ji = 1, jpi 
    196                         zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 
    197                         irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
    198                         zekb(ji,jj) = rkrgb(1,irgb) 
    199                         zekg(ji,jj) = rkrgb(2,irgb) 
    200                         zekr(ji,jj) = rkrgb(3,irgb) 
    201                      END DO 
    202                   END DO 
    203                ELSE                                            ! Variable ocean volume but constant chrlorophyll 
    204                   zchl = 0.05                                     ! constant chlorophyll 
    205                   irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 
    206                   zekb(:,:) = rkrgb(1,irgb)                       ! Separation in R-G-B depending of the chlorophyll  
    207                   zekg(:,:) = rkrgb(2,irgb) 
    208                   zekr(:,:) = rkrgb(3,irgb) 
     203                     DO ji = 1, jpi 
     204                        zchl    = sf_chl(1)%fnow(ji,jj,1) 
     205                        zCtot   = 40.6  * zchl**0.459 
     206                        zze     = 568.2 * zCtot**(-0.746) 
     207                        IF( zze > 102. ) zze = 200.0 * zCtot**(-0.293) 
     208                        zlogc   = LOG( zchl ) 
     209                        zlogc2  = zlogc * zlogc 
     210                        zlogc3  = zlogc * zlogc * zlogc 
     211                        zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3 
     212                        zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2 
     213                        zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 
     214                        zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 
     215                        zCze    = 1.12  * (zchl)**0.803  
     216                        DO jk = 1, nksr + 1 
     217                           zpsi = fsdept(ji,jj,jk) / zze 
     218                           zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) ) 
     219                        END DO 
     220                        ! 
     221                      END DO 
     222                   END DO 
     223                     ! 
     224               ELSE                              !* Variable ocean volume but constant chrlorophyll 
     225                  DO jk = 1, nksr + 1 
     226                     zchl3d(:,:,jk) = 0.05  
     227                  ENDDO 
    209228               ENDIF 
    210229               ! 
    211                zcoef  = ( 1. - rn_abs ) / 3.e0                        ! equi-partition in R-G-B 
     230               zcoef  = ( 1. - rn_abs ) / 3.e0                        !  equi-partition in R-G-B 
    212231               ze0(:,:,1) = rn_abs  * qsr(:,:) 
    213232               ze1(:,:,1) = zcoef * qsr(:,:) 
     
    217236               ! 
    218237               DO jk = 2, nksr+1 
     238                  ! 
     239                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of vertical profile of Chl 
     240!CDIR NOVERRCHK 
     241                     DO ji = 1, jpi 
     242                        zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,jk) ) ) 
     243                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
     244                        zekb(ji,jj) = rkrgb(1,irgb) 
     245                        zekg(ji,jj) = rkrgb(2,irgb) 
     246                        zekr(ji,jj) = rkrgb(3,irgb) 
     247                     END DO 
     248                  END DO 
    219249!CDIR NOVERRCHK 
    220250                  DO jj = 1, jpj 
     
    233263                  END DO 
    234264               END DO 
    235                ! clem: store attenuation coefficient of the first ocean level 
    236                IF ( ln_qsr_ice ) THEN 
    237                   DO jj = 1, jpj 
    238                      DO ji = 1, jpi 
    239                         zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     ) 
    240                         zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 
    241                         zzc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 
    242                         zzc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 
    243                         fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 + zzc2  + zzc3  ) * tmask(ji,jj,2)  
    244                      END DO 
    245                   END DO 
    246                ENDIF 
    247265               ! 
    248266               DO jk = 1, nksr                                        ! compute and add qsr trend to ta 
     
    251269               zea(:,:,nksr+1:jpk) = 0.e0     ! below 400m set to zero 
    252270               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution 
     271               ! 
     272               IF ( ln_qsr_ice ) THEN    ! store attenuation coefficient of the first ocean level 
     273!CDIR NOVERRCHK 
     274                  DO jj = 1, jpj                                         ! Separation in R-G-B depending of the surface Chl 
     275!CDIR NOVERRCHK 
     276                     DO ji = 1, jpi 
     277                        zchl = MIN( 10. , MAX( 0.03, zchl3d(ji,jj,1) ) ) 
     278                        irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
     279                        zekb(ji,jj) = rkrgb(1,irgb) 
     280                        zekg(ji,jj) = rkrgb(2,irgb) 
     281                        zekr(ji,jj) = rkrgb(3,irgb) 
     282                     END DO 
     283                  END DO 
     284                  !  
     285                  DO jj = 1, jpj 
     286                     DO ji = 1, jpi 
     287                        zc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     ) 
     288                        zc1 = zcoef  * EXP( - fse3t(ji,jj,1) * zekb(ji,jj) ) 
     289                        zc2 = zcoef  * EXP( - fse3t(ji,jj,1) * zekg(ji,jj) ) 
     290                        zc3 = zcoef  * EXP( - fse3t(ji,jj,1) * zekr(ji,jj) ) 
     291                        fraqsr_1lev(ji,jj) = 1.0 - ( zc0 + zc1 + zc2  + zc3  ) * tmask(ji,jj,2)  
     292                     END DO 
     293                  END DO 
     294                  ! 
     295               ENDIF 
    253296               ! 
    254297            ELSE                                                 !*  Constant Chlorophyll 
     
    256299                  qsr_hc(:,:,jk) =  etot3(:,:,jk) * qsr(:,:) 
    257300               END DO 
    258                ! clem: store attenuation coefficient of the first ocean level 
    259                IF ( ln_qsr_ice ) THEN 
     301               ! store attenuation coefficient of the first ocean level 
     302               IF( ln_qsr_ice ) THEN 
    260303                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    261304               ENDIF 
     
    339382      ! 
    340383      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    341       CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  
     384      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea, zchl3d )  
    342385      ! 
    343386      IF( nn_timing == 1 )  CALL timing_stop('tra_qsr') 
     
    405448         WRITE(numout,*) '      bio-model            light penetration   ln_qsr_bio = ', ln_qsr_bio 
    406449         WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice 
    407          WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)    nn_chldta  = ', nn_chldta 
     450         WRITE(numout,*) '      RGB : Chl data (=1/2) or cst value (=0)  nn_chldta  = ', nn_chldta 
    408451         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs = ', rn_abs 
    409452         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0 
     
    429472         IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr =  1  
    430473         IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr =  2 
    431          IF( ln_qsr_2bd                      )   nqsr =  3 
    432          IF( ln_qsr_bio                      )   nqsr =  4 
     474         IF( ln_qsr_rgb .AND. nn_chldta == 2 )   nqsr =  3 
     475         IF( ln_qsr_2bd                      )   nqsr =  4 
     476         IF( ln_qsr_bio                      )   nqsr =  5 
    433477         ! 
    434478         IF(lwp) THEN                   ! Print the choice 
    435479            WRITE(numout,*) 
    436480            IF( nqsr ==  1 )   WRITE(numout,*) '         R-G-B   light penetration - Constant Chlorophyll' 
    437             IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - Chl data ' 
    438             IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration' 
    439             IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration' 
     481            IF( nqsr ==  2 )   WRITE(numout,*) '         R-G-B   light penetration - 2D Chl data ' 
     482            IF( nqsr ==  3 )   WRITE(numout,*) '         R-G-B   light penetration - 3D Chl data ' 
     483            IF( nqsr ==  4 )   WRITE(numout,*) '         2 bands light penetration' 
     484            IF( nqsr ==  5 )   WRITE(numout,*) '         bio-model light penetration' 
    440485         ENDIF 
    441486         ! 
     
    460505            IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
    461506            ! 
    462             IF( nn_chldta == 1 ) THEN           !* Chl data : set sf_chl structure 
     507            IF( nn_chldta == 1  .OR. nn_chldta == 2 ) THEN           !* Chl data : set sf_chl structure 
    463508               IF(lwp) WRITE(numout,*) 
    464509               IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file' 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r5781 r6333  
    177177                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  ) 
    178178               ! add to the eddy viscosity coef. previously computed 
     179# if defined key_zdftmx_new 
     180               ! key_zdftmx_new: New internal wave-driven param: use avs value computed by zdftmx 
     181               avs (ji,jj,jk) = avs(ji,jj,jk) + zavfs + zavds 
     182# else 
    179183               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds 
     184# endif 
    180185               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt 
    181186               avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds ) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r6237 r6333  
    357357            DO ji = fs_2, fs_jpim1   ! vector opt. 
    358358               zcof   = zfact1 * tmask(ji,jj,jk) 
     359# if defined key_zdftmx_new 
     360               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability) 
     361               zzd_up = zcof * ( MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) )   &  ! upper diagonal 
     362                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) ) 
     363               zzd_lw = zcof * ( MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) )   &  ! lower diagonal 
     364                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) ) 
     365# else 
    359366               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal 
    360367                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) ) 
    361368               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal 
    362369                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) ) 
     370# endif 
    363371                  !                                                           ! shear prod. at w-point weightened by mask 
    364372               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
     
    735743      ! 
    736744      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number 
     745# if defined key_zdftmx_new 
     746      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used 
     747      rn_emin  = 1.e-10_wp 
     748      rmxl_min = 1.e-03_wp 
     749      IF(lwp) THEN                  ! Control print 
     750         WRITE(numout,*) 
     751         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 ' 
     752         WRITE(numout,*) '~~~~~~~~~~~~' 
     753      ENDIF 
     754# else 
    737755      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity 
     756# endif 
    738757      ! 
    739758      IF(lwp) THEN                    !* Control print 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r5781 r6333  
    561561   END SUBROUTINE zdf_tmx_init 
    562562 
     563#elif defined key_zdftmx_new 
     564   !!---------------------------------------------------------------------- 
     565   !!   'key_zdftmx_new'               Internal wave-driven vertical mixing 
     566   !!---------------------------------------------------------------------- 
     567   !!   zdf_tmx       : global     momentum & tracer Kz with wave induced Kz 
     568   !!   zdf_tmx_init  : global     momentum & tracer Kz with wave induced Kz 
     569   !!---------------------------------------------------------------------- 
     570   USE oce            ! ocean dynamics and tracers variables 
     571   USE dom_oce        ! ocean space and time domain variables 
     572   USE zdf_oce        ! ocean vertical physics variables 
     573   USE zdfddm         ! ocean vertical physics: double diffusive mixing 
     574   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     575   USE eosbn2         ! ocean equation of state 
     576   USE phycst         ! physical constants 
     577   USE prtctl         ! Print control 
     578   USE in_out_manager ! I/O manager 
     579   USE iom            ! I/O Manager 
     580   USE lib_mpp        ! MPP library 
     581   USE wrk_nemo       ! work arrays 
     582   USE timing         ! Timing 
     583   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     584 
     585   IMPLICIT NONE 
     586   PRIVATE 
     587 
     588   PUBLIC   zdf_tmx         ! called in step module  
     589   PUBLIC   zdf_tmx_init    ! called in nemogcm module  
     590   PUBLIC   zdf_tmx_alloc   ! called in nemogcm module 
     591 
     592   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: wave-driven mixing flag 
     593 
     594   !                       !!* Namelist  namzdf_tmx : internal wave-driven mixing * 
     595   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2) 
     596   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency 
     597   LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F) 
     598 
     599   REAL(wp) ::  r1_6 = 1._wp / 6._wp 
     600 
     601   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_tmx     ! power available from high-mode wave breaking (W/m2) 
     602   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_tmx     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2) 
     603   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_tmx     ! power available from low-mode, critical slope wave breaking (W/m2) 
     604   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_tmx     ! WKB decay scale for high-mode energy dissipation (m) 
     605   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_tmx     ! decay scale for low-mode critical slope dissipation (m) 
     606   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_tmx     ! local energy density available for mixing (W/kg) 
     607   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_tmx     ! buoyancy flux Kz * N^2 (W/kg) 
     608   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_tmx    ! vertically integrated buoyancy flux (W/m2) 
     609   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T) 
     610   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity 
     611 
     612   !! * Substitutions 
     613#  include "zdfddm_substitute.h90" 
     614#  include "domzgr_substitute.h90" 
     615#  include "vectopt_loop_substitute.h90" 
     616   !!---------------------------------------------------------------------- 
     617   !! NEMO/OPA 4.0 , NEMO Consortium (2016) 
     618   !! $Id$ 
     619   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     620   !!---------------------------------------------------------------------- 
     621CONTAINS 
     622 
     623   INTEGER FUNCTION zdf_tmx_alloc() 
     624      !!---------------------------------------------------------------------- 
     625      !!                ***  FUNCTION zdf_tmx_alloc  *** 
     626      !!---------------------------------------------------------------------- 
     627      ALLOCATE(     ebot_tmx(jpi,jpj),  epyc_tmx(jpi,jpj),  ecri_tmx(jpi,jpj)    ,   & 
     628      &             hbot_tmx(jpi,jpj),  hcri_tmx(jpi,jpj),  emix_tmx(jpi,jpj,jpk),   & 
     629      &         bflx_tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk),   &  
     630      &         zav_wave(jpi,jpj,jpk), STAT=zdf_tmx_alloc     ) 
     631      ! 
     632      IF( lk_mpp             )   CALL mpp_sum ( zdf_tmx_alloc ) 
     633      IF( zdf_tmx_alloc /= 0 )   CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays') 
     634   END FUNCTION zdf_tmx_alloc 
     635 
     636 
     637   SUBROUTINE zdf_tmx( kt ) 
     638      !!---------------------------------------------------------------------- 
     639      !!                  ***  ROUTINE zdf_tmx  *** 
     640      !!                    
     641      !! ** Purpose :   add to the vertical mixing coefficients the effect of 
     642      !!              breaking internal waves. 
     643      !! 
     644      !! ** Method  : - internal wave-driven vertical mixing is given by: 
     645      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_tmx /( Nu * N^2 )  ) 
     646      !!              where emix_tmx is the 3D space distribution of the wave-breaking  
     647      !!              energy and Nu the molecular kinematic viscosity. 
     648      !!              The function f(Reb) is linear (constant mixing efficiency) 
     649      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T. 
     650      !! 
     651      !!              - Compute emix_tmx, the 3D power density that allows to compute 
     652      !!              Reb and therefrom the wave-induced vertical diffusivity. 
     653      !!              This is divided into three components: 
     654      !!                 1. Bottom-intensified low-mode dissipation at critical slopes 
     655      !!                     emix_tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx ) 
     656      !!                                   / ( 1. - EXP( - H/hcri_tmx ) ) * hcri_tmx 
     657      !!              where hcri_tmx is the characteristic length scale of the bottom  
     658      !!              intensification, ecri_tmx a map of available power, and H the ocean depth. 
     659      !!                 2. Pycnocline-intensified low-mode dissipation 
     660      !!                     emix_tmx(z) = ( epyc_tmx / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
     661      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
     662      !!              where epyc_tmx is a map of available power, and nn_zpyc 
     663      !!              is the chosen stratification-dependence of the internal wave 
     664      !!              energy dissipation. 
     665      !!                 3. WKB-height dependent high mode dissipation 
     666      !!                     emix_tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx) 
     667      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_tmx) * e3w(z) ) 
     668      !!              where hbot_tmx is the characteristic length scale of the WKB bottom  
     669      !!              intensification, ebot_tmx is a map of available power, and z_wkb is the 
     670      !!              WKB-stretched height above bottom defined as 
     671      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) ) 
     672      !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    ) 
     673      !! 
     674      !!              - update the model vertical eddy viscosity and diffusivity:  
     675      !!                     avt  = avt  +    av_wave 
     676      !!                     avm  = avm  +    av_wave 
     677      !!                     avmu = avmu + mi(av_wave) 
     678      !!                     avmv = avmv + mj(av_wave) 
     679      !! 
     680      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing: 
     681      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb) 
     682      !! 
     683      !! ** Action  : - Define emix_tmx used to compute internal wave-induced mixing 
     684      !!              - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing     
     685      !! 
     686      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep. 
     687      !!---------------------------------------------------------------------- 
     688      INTEGER, INTENT(in) ::   kt   ! ocean time-step  
     689      ! 
     690      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     691      REAL(wp) ::   ztpc         ! scalar workspace 
     692      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure 
     693      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth 
     694      REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom 
     695      REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution 
     696      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid) 
     697      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid) 
     698      REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter 
     699      !!---------------------------------------------------------------------- 
     700      ! 
     701      IF( nn_timing == 1 )   CALL timing_start('zdf_tmx') 
     702      ! 
     703      CALL wrk_alloc( jpi,jpj,       zfact, zhdep ) 
     704      CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
     705 
     706      !                          ! ----------------------------- ! 
     707      !                          !  Internal wave-driven mixing  !  (compute zav_wave) 
     708      !                          ! ----------------------------- ! 
     709      !                              
     710      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
     711      !                                                 using an exponential decay from the seafloor. 
     712      DO jj = 1, jpj                ! part independent of the level 
     713         DO ji = 1, jpi 
     714            zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
     715            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_tmx(ji,jj) )  ) 
     716            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj) 
     717         END DO 
     718      END DO 
     719 
     720      DO jk = 2, jpkm1              ! complete with the level-dependent part 
     721         emix_tmx(:,:,jk) = zfact(:,:) * (  EXP( ( fsde3w(:,:,jk  ) - zhdep(:,:) ) / hcri_tmx(:,:) )                      & 
     722            &                             - EXP( ( fsde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_tmx(:,:) )  ) * wmask(:,:,jk)   & 
     723            &                          / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) ) 
     724      END DO 
     725 
     726      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying  
     727      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 
     728 
     729      SELECT CASE ( nn_zpyc ) 
     730 
     731      CASE ( 1 )               ! Dissipation scales as N (recommended) 
     732 
     733         zfact(:,:) = 0._wp 
     734         DO jk = 2, jpkm1              ! part independent of the level 
     735            zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     736         END DO 
     737 
     738         DO jj = 1, jpj 
     739            DO ji = 1, jpi 
     740               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     741            END DO 
     742         END DO 
     743 
     744         DO jk = 2, jpkm1              ! complete with the level-dependent part 
     745            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     746         END DO 
     747 
     748      CASE ( 2 )               ! Dissipation scales as N^2 
     749 
     750         zfact(:,:) = 0._wp 
     751         DO jk = 2, jpkm1              ! part independent of the level 
     752            zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     753         END DO 
     754 
     755         DO jj= 1, jpj 
     756            DO ji = 1, jpi 
     757               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     758            END DO 
     759         END DO 
     760 
     761         DO jk = 2, jpkm1              ! complete with the level-dependent part 
     762            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     763         END DO 
     764 
     765      END SELECT 
     766 
     767      !                        !* WKB-height dependent mixing: distribute energy over the time-varying  
     768      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
     769       
     770      zwkb(:,:,:) = 0._wp 
     771      zfact(:,:) = 0._wp 
     772      DO jk = 2, jpkm1 
     773         zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     774         zwkb(:,:,jk) = zfact(:,:) 
     775      END DO 
     776 
     777      DO jk = 2, jpkm1 
     778         DO jj = 1, jpj 
     779            DO ji = 1, jpi 
     780               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
     781                                            &           * tmask(ji,jj,jk) / zfact(ji,jj) 
     782            END DO 
     783         END DO 
     784      END DO 
     785      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 
     786 
     787      zweight(:,:,:) = 0._wp 
     788      DO jk = 2, jpkm1 
     789         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_tmx(:,:) * wmask(:,:,jk)                    & 
     790            &   * (  EXP( -zwkb(:,:,jk) / hbot_tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) )  ) 
     791      END DO 
     792 
     793      zfact(:,:) = 0._wp 
     794      DO jk = 2, jpkm1              ! part independent of the level 
     795         zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
     796      END DO 
     797 
     798      DO jj = 1, jpj 
     799         DO ji = 1, jpi 
     800            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     801         END DO 
     802      END DO 
     803 
     804      DO jk = 2, jpkm1              ! complete with the level-dependent part 
     805         emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
     806            &                                / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) ) 
     807      END DO 
     808 
     809 
     810      ! Calculate molecular kinematic viscosity 
     811      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
     812         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0 
     813      DO jk = 2, jpkm1 
     814         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
     815      END DO 
     816 
     817      ! Calculate turbulence intensity parameter Reb 
     818      DO jk = 2, jpkm1 
     819         zReb(:,:,jk) = emix_tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
     820      END DO 
     821 
     822      ! Define internal wave-induced diffusivity 
     823      DO jk = 2, jpkm1 
     824         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
     825      END DO 
     826 
     827      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
     828         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
     829            DO jj = 1, jpj 
     830               DO ji = 1, jpi 
     831                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
     832                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     833                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 
     834                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     835                  ENDIF 
     836               END DO 
     837            END DO 
     838         END DO 
     839      ENDIF 
     840 
     841      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
     842         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
     843      END DO 
     844 
     845      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
     846         ztpc = 0._wp 
     847         DO jk = 2, jpkm1 
     848            DO jj = 1, jpj 
     849               DO ji = 1, jpi 
     850                  ztpc = ztpc + fse3w(ji,jj,jk) * e1e2t(ji,jj)   & 
     851                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
     852               END DO 
     853            END DO 
     854         END DO 
     855         IF( lk_mpp )   CALL mpp_sum( ztpc ) 
     856         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
     857  
     858         IF(lwp) THEN 
     859            WRITE(numout,*) 
     860            WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)' 
     861            WRITE(numout,*) '~~~~~~~ ' 
     862            WRITE(numout,*) 
     863            WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW' 
     864         ENDIF 
     865      ENDIF 
     866 
     867      !                          ! ----------------------- ! 
     868      !                          !   Update  mixing coefs  !                           
     869      !                          ! ----------------------- ! 
     870      !       
     871      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
     872         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb 
     873            DO jj = 1, jpj 
     874               DO ji = 1, jpi 
     875                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  & 
     876                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   & 
     877                      &                 ) * wmask(ji,jj,jk) 
     878               END DO 
     879            END DO 
     880         END DO 
     881         CALL iom_put( "av_ratio", zav_ratio ) 
     882         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
     883            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
     884            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     885            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     886         END DO 
     887         ! 
     888      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
     889         DO jk = 2, jpkm1 
     890            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     891            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     892            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     893         END DO 
     894      ENDIF 
     895 
     896      DO jk = 2, jpkm1              !* update momentum diffusivity at wu and wv points 
     897         DO jj = 2, jpjm1 
     898            DO ji = fs_2, fs_jpim1  ! vector opt. 
     899               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     900               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
     901            END DO 
     902         END DO 
     903      END DO 
     904      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition 
     905 
     906      !                             !* output internal wave-driven mixing coefficient 
     907      CALL iom_put( "av_wave", zav_wave ) 
     908                                    !* output useful diagnostics: N^2, Kz * N^2 (bflx_tmx),  
     909                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_tmx), energy density (emix_tmx) 
     910      IF( iom_use("bflx_tmx") .OR. iom_use("pcmap_tmx") ) THEN 
     911         bflx_tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
     912         pcmap_tmx(:,:) = 0._wp 
     913         DO jk = 2, jpkm1 
     914            pcmap_tmx(:,:) = pcmap_tmx(:,:) + fse3w(:,:,jk) * bflx_tmx(:,:,jk) * wmask(:,:,jk) 
     915         END DO 
     916         pcmap_tmx(:,:) = rau0 * pcmap_tmx(:,:) 
     917         CALL iom_put( "bflx_tmx", bflx_tmx ) 
     918         CALL iom_put( "pcmap_tmx", pcmap_tmx ) 
     919      ENDIF 
     920      CALL iom_put( "bn2", rn2 ) 
     921      CALL iom_put( "emix_tmx", emix_tmx ) 
     922       
     923      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep ) 
     924      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
     925 
     926      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 
     927      ! 
     928      IF( nn_timing == 1 )   CALL timing_stop('zdf_tmx') 
     929      ! 
     930   END SUBROUTINE zdf_tmx 
     931 
     932 
     933   SUBROUTINE zdf_tmx_init 
     934      !!---------------------------------------------------------------------- 
     935      !!                  ***  ROUTINE zdf_tmx_init  *** 
     936      !!                      
     937      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading 
     938      !!              of input power maps and decay length scales in netcdf files. 
     939      !! 
     940      !! ** Method  : - Read the namzdf_tmx namelist and check the parameters 
     941      !! 
     942      !!              - Read the input data in NetCDF files : 
     943      !!              power available from high-mode wave breaking (mixing_power_bot.nc) 
     944      !!              power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc) 
     945      !!              power available from critical slope wave-breaking (mixing_power_cri.nc) 
     946      !!              WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc) 
     947      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc) 
     948      !! 
     949      !! ** input   : - Namlist namzdf_tmx 
     950      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc, 
     951      !!              decay_scale_bot.nc decay_scale_cri.nc 
     952      !! 
     953      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter 
     954      !!              - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx 
     955      !! 
     956      !! References : de Lavergne et al. 2015, JPO; 2016, in prep. 
     957      !!          
     958      !!---------------------------------------------------------------------- 
     959      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     960      INTEGER  ::   inum         ! local integer 
     961      INTEGER  ::   ios 
     962      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars 
     963      !! 
     964      NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff 
     965      !!---------------------------------------------------------------------- 
     966      ! 
     967      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init') 
     968      ! 
     969      REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing 
     970      READ  ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901) 
     971901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp ) 
     972      ! 
     973      REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing 
     974      READ  ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 ) 
     975902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp ) 
     976      IF(lwm) WRITE ( numond, namzdf_tmx_new ) 
     977      ! 
     978      IF(lwp) THEN                  ! Control print 
     979         WRITE(numout,*) 
     980         WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing' 
     981         WRITE(numout,*) '~~~~~~~~~~~~' 
     982         WRITE(numout,*) '   Namelist namzdf_tmx_new : set wave-driven mixing parameters' 
     983         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc 
     984         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar 
     985         WRITE(numout,*) '      Differential internal wave-driven mixing (T) or not (F)   = ', ln_tsdiff 
     986      ENDIF 
     987       
     988      ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and 
     989      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should  
     990      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6). 
     991      avmb(:) = 1.4e-6_wp        ! viscous molecular value 
     992      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_tmx)     
     993      avtb_2d(:,:) = 1.e0_wp     ! uniform  
     994      IF(lwp) THEN                  ! Control print 
     995         WRITE(numout,*) 
     996         WRITE(numout,*) '   Force the background value applied to avm & avt in TKE to be everywhere ',   & 
     997            &               'the viscous molecular value & a very small diffusive value, resp.' 
     998      ENDIF 
     999       
     1000      IF( .NOT.lk_zdfddm )   CALL ctl_stop( 'STOP', 'zdf_tmx_init_new : key_zdftmx_new requires key_zdfddm' ) 
     1001       
     1002      !                             ! allocate tmx arrays 
     1003      IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' ) 
     1004      ! 
     1005      !                             ! read necessary fields 
     1006      CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2] 
     1007      CALL iom_get  (inum, jpdom_data, 'field', ebot_tmx, 1 )  
     1008      CALL iom_close(inum) 
     1009      ! 
     1010      CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2] 
     1011      CALL iom_get  (inum, jpdom_data, 'field', epyc_tmx, 1 ) 
     1012      CALL iom_close(inum) 
     1013      ! 
     1014      CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2] 
     1015      CALL iom_get  (inum, jpdom_data, 'field', ecri_tmx, 1 ) 
     1016      CALL iom_close(inum) 
     1017      ! 
     1018      CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m] 
     1019      CALL iom_get  (inum, jpdom_data, 'field', hbot_tmx, 1 ) 
     1020      CALL iom_close(inum) 
     1021      ! 
     1022      CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m] 
     1023      CALL iom_get  (inum, jpdom_data, 'field', hcri_tmx, 1 ) 
     1024      CALL iom_close(inum) 
     1025 
     1026      ebot_tmx(:,:) = ebot_tmx(:,:) * ssmask(:,:) 
     1027      epyc_tmx(:,:) = epyc_tmx(:,:) * ssmask(:,:) 
     1028      ecri_tmx(:,:) = ecri_tmx(:,:) * ssmask(:,:) 
     1029 
     1030      ! Set once for all to zero the first and last vertical levels of appropriate variables 
     1031      emix_tmx (:,:, 1 ) = 0._wp 
     1032      emix_tmx (:,:,jpk) = 0._wp 
     1033      zav_ratio(:,:, 1 ) = 0._wp 
     1034      zav_ratio(:,:,jpk) = 0._wp 
     1035      zav_wave (:,:, 1 ) = 0._wp 
     1036      zav_wave (:,:,jpk) = 0._wp 
     1037 
     1038      zbot = glob_sum( e1e2t(:,:) * ebot_tmx(:,:) ) 
     1039      zpyc = glob_sum( e1e2t(:,:) * epyc_tmx(:,:) ) 
     1040      zcri = glob_sum( e1e2t(:,:) * ecri_tmx(:,:) ) 
     1041      IF(lwp) THEN 
     1042         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW' 
     1043         WRITE(numout,*) '      Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW' 
     1044         WRITE(numout,*) '      Critical slope wave-breaking energy:        ', zcri * 1.e-12_wp, 'TW' 
     1045      ENDIF 
     1046      ! 
     1047      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init') 
     1048      ! 
     1049   END SUBROUTINE zdf_tmx_init 
     1050 
    5631051#else 
    5641052   !!---------------------------------------------------------------------- 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90

    r6237 r6333  
    3131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sio3eq   ! chemistry of Si 
    3232   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fekeq    ! chemistry of Fe 
    33    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemc    ! Solubilities of O2 and CO2 
    34    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2    ! Solubilities of O2 and CO2 
     33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ::   chemc    ! Solubilities of O2 and CO2 
     34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   chemo2   ! Solubilities of O2 and CO2 
    3535 
    3636   REAL(wp), PUBLIC ::   atcox  = 0.20946         ! units atm 
     
    7676   REAL(wp) ::   st1    =      0.14     ! constants for calculate concentrations for sulfate 
    7777   REAL(wp) ::   st2    =  1./96.062    !  (Morris & Riley 1966) 
    78    REAL(wp) ::   ks0    =    141.328  
    79    REAL(wp) ::   ks1    =  -4276.1   
    80    REAL(wp) ::   ks2    =    -23.093 
    81    REAL(wp) ::   ks3    = -13856.   
    82    REAL(wp) ::   ks4    =   324.57  
    83    REAL(wp) ::   ks5    =   -47.986 
    84    REAL(wp) ::   ks6    =  35474.  
    85    REAL(wp) ::   ks7    =   -771.54 
    86    REAL(wp) ::   ks8    =    114.723 
    87    REAL(wp) ::   ks9    =  -2698.   
    88    REAL(wp) ::   ks10   =   1776.  
    89    REAL(wp) ::   ks11   =      1. 
    90    REAL(wp) ::   ks12   =     -0.001005  
    9178 
    9279   REAL(wp) ::   ft1    =    0.000067   ! constants for calculate concentrations for fluorides 
    9380   REAL(wp) ::   ft2    = 1./18.9984    ! (Dickson & Riley 1979 ) 
    94    REAL(wp) ::   kf0    =  -12.641     
    95    REAL(wp) ::   kf1    = 1590.2     
    96    REAL(wp) ::   kf2    =    1.525     
    97    REAL(wp) ::   kf3    =    1.0      
    98    REAL(wp) ::   kf4    =   -0.001005 
    99  
    100    REAL(wp) ::   cb0    = -8966.90      ! Coeff. for 1. dissoc. of boric acid  
    101    REAL(wp) ::   cb1    = -2890.53      ! (Dickson and Goyet, 1994) 
    102    REAL(wp) ::   cb2    =   -77.942 
    103    REAL(wp) ::   cb3    =     1.728 
    104    REAL(wp) ::   cb4    =    -0.0996 
    105    REAL(wp) ::   cb5    =   148.0248 
    106    REAL(wp) ::   cb6    =   137.1942 
    107    REAL(wp) ::   cb7    =     1.62142 
    108    REAL(wp) ::   cb8    =   -24.4344 
    109    REAL(wp) ::   cb9    =   -25.085 
    110    REAL(wp) ::   cb10   =    -0.2474  
    111    REAL(wp) ::   cb11   =     0.053105 
    112  
    113    REAL(wp) ::   cw0    = -13847.26     ! Coeff. for dissoc. of water (Dickson and Riley, 1979 ) 
    114    REAL(wp) ::   cw1    =    148.9652   
    115    REAL(wp) ::   cw2    =    -23.6521 
    116    REAL(wp) ::   cw3    =    118.67  
    117    REAL(wp) ::   cw4    =     -5.977  
    118    REAL(wp) ::   cw5    =      1.0495   
    119    REAL(wp) ::   cw6    =     -0.01615 
    12081 
    12182   !                                    ! volumetric solubility constants for o2 in ml/L   
     
    200161         DO ji = 1, jpi 
    201162            !                             ! SET ABSOLUTE TEMPERATURE 
    202             ztkel = tsn(ji,jj,1,jp_tem) + 273.16 
     163            ztkel = tsn(ji,jj,1,jp_tem) + 273.15 
    203164            zt    = ztkel * 0.01 
    204165            zt2   = zt * zt 
     
    209170            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 
    210171            zcek1 = ca0 + ca1 / zt + ca2 * zlogt + ca3 * zt2 + zsal * ( ca4 + ca5 * zt + ca6 * zt2 ) 
    211             !                             ! LN(K0) OF SOLUBILITY OF O2 and N2 in ml/L (EQ. 8, GARCIA AND GORDON, 1992) 
    212             ztgg  = LOG( ( 298.15 - tsn(ji,jj,1,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature 
    213             ztgg2 = ztgg  * ztgg 
    214             ztgg3 = ztgg2 * ztgg 
    215             ztgg4 = ztgg3 * ztgg 
    216             ztgg5 = ztgg4 * ztgg 
    217             zoxy  = ox0 + ox1 * ztgg + ox2 * ztgg2 + ox3 * ztgg3 + ox4 * ztgg4 + ox5 * ztgg5   & 
    218                    + zsal * ( ox6 + ox7 * ztgg + ox8 * ztgg2 + ox9 * ztgg3 ) +  ox10 * zsal2 
    219  
    220172            !                             ! SET SOLUBILITIES OF O2 AND CO2  
    221             chemc(ji,jj,1) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm) 
    222             chemc(ji,jj,2) = ( EXP( zoxy  ) * o2atm ) * oxyco              ! mol/(L atm) 
     173            chemc(ji,jj) = EXP( zcek1 ) * 1.e-6 * rhop(ji,jj,1) / 1000.  ! mol/(L uatm) 
    223174            ! 
    224175         END DO 
     
    233184!CDIR NOVERRCHK 
    234185            DO ji = 1, jpi 
    235               ztkel = tsn(ji,jj,jk,jp_tem) + 273.16 
     186              ztkel = tsn(ji,jj,jk,jp_tem) + 273.15 
    236187              zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35. 
    237188              zsal2 = zsal * zsal 
     
    263214 
    264215               ! SET ABSOLUTE TEMPERATURE 
    265                ztkel   = tsn(ji,jj,jk,jp_tem) + 273.16 
     216               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.15 
    266217               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35. 
    267218               zsqrt  = SQRT( zsal ) 
     
    284235 
    285236               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 
    286                zcks    = EXP(  ks1 * ztr + ks0 + ks2 * zlogt                           & 
    287                   &                     + ( ks3 * ztr + ks4 + ks5 * zlogt ) * zisqrt   & 
    288                   &                     + ( ks6 * ztr + ks7 + ks8 * zlogt ) * zis      & 
    289                   &                     + ks9 * ztr * zis * zisqrt + ks10 * ztr *zis2 + LOG( ks11 + ks12 *zsal )  ) 
     237               zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         & 
     238               &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt & 
     239               &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    & 
     240               &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         & 
     241               &         + LOG(1.0 - 0.001005 * zsal)) 
     242               ! 
     243               aphscale(ji,jj,jk) = ( 1. + zst / zcks ) 
    290244 
    291245               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 
    292                zckf    = EXP(  kf1 * ztr + kf0 + kf2 * zisqrt + LOG( kf3 + kf4 * zsal )  ) 
     246               zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   & 
     247               &         + LOG(1.0d0 - 0.001005d0*zsal)            & 
     248               &         + LOG(1.0d0 + zst/zcks)) 
    293249 
    294250               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 
    295                zckb    = ( cb0 + cb1 * zsqrt + cb2  * zsal + cb3 * zsal15 + cb4 * zsal * zsal ) * ztr   & 
    296                   &    + ( cb5 + cb6 * zsqrt + cb7  * zsal )                                            & 
    297                   &    + ( cb8 + cb9 * zsqrt + cb10 * zsal ) * zlogt + cb11 * zsqrt * ztkel             & 
    298                   &    + LOG(  ( 1.+ zst / zcks + zft / zckf ) / ( 1.+ zst / zcks )  ) 
     251               zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        & 
     252               &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         & 
     253               &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   & 
     254               &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      &  
     255               &      * zlogt + 0.053105*zsqrt*ztkel 
     256 
    299257 
    300258               zck1    = c10 * ztr + c11 + c12 * zlogt + c13 * zsal + c14 * zsal * zsal 
     
    302260 
    303261               ! PKW (H2O) (DICKSON AND RILEY, 1979) 
    304                zckw    = cw0 * ztr + cw1 + cw2 * zlogt + ( cw3 * ztr + cw4 + cw5 * zlogt ) * zsqrt + cw6 * zsal 
    305  
     262               zckw = -13847.26*ztr + 148.9652 - 23.6521 * zlogt    &  
     263               &     + (118.67*ztr - 5.977 + 1.0495 * zlogt)        & 
     264               &     * zsqrt - 0.01615 * zsal 
    306265 
    307266               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 
     
    378337      !!                     ***  ROUTINE p4z_che_alloc  *** 
    379338      !!---------------------------------------------------------------------- 
    380       ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,2), chemo2(jpi,jpj,jpk), STAT=p4z_che_alloc ) 
     339      ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj), chemo2(jpi,jpj,jpk),   & 
     340      &         STAT=p4z_che_alloc ) 
    381341      ! 
    382342      IF( p4z_che_alloc /= 0 )   CALL ctl_warn('p4z_che_alloc : failed to allocate arrays.') 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90

    r6237 r6333  
    8484      ! 
    8585      INTEGER  ::   ji, jj, jm, iind, iindm1 
    86       REAL(wp) ::   ztc, ztc2, ztc3, zws, zkgwan 
     86      REAL(wp) ::   ztc, ztc2, ztc3, ztc4, zws, zkgwan 
    8787      REAL(wp) ::   zfld, zflu, zfld16, zflu16, zfact 
    8888      REAL(wp) ::   zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 
     
    135135 
    136136               ! CALCULATE [ALK]([CO3--], [HCO3-]) 
    137                zalk  = zalka - (  akw3(ji,jj,1) / zph - zph + zbot / ( 1.+ zph / akb3(ji,jj,1) )  ) 
     137               zalk  = zalka - (  akw3(ji,jj,1) / zph - zph / aphscale(ji,jj,1)    & 
     138               &       + zbot / ( 1.+ zph / akb3(ji,jj,1) )  ) 
    138139 
    139140               ! CALCULATE [H+] AND [H2CO3] 
     
    162163            ztc2 = ztc * ztc 
    163164            ztc3 = ztc * ztc2  
     165            ztc4 = ztc2 * ztc2  
    164166            ! Compute the schmidt Number both O2 and CO2 
    165             zsch_co2 = 2073.1 - 125.62 * ztc + 3.6276 * ztc2 - 0.043126 * ztc3 
    166             zsch_o2  = 1953.4 - 128.0  * ztc + 3.9918 * ztc2 - 0.050091 * ztc3 
     167            zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4 
     168            zsch_o2  = 1920.4 - 135.6  * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4 
    167169            !  wind speed  
    168170            zws  = wndm(ji,jj) * wndm(ji,jj) 
    169171            ! Compute the piston velocity for O2 and CO2 
    170             zkgwan = 0.3 * zws  + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946  * ztc2 ) 
     172            zkgwan = 0.251 * zws 
    171173            zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 
    172174# if defined key_degrad 
     
    182184         DO ji = 1, jpi 
    183185            ! Compute CO2 flux for the sea and air 
    184             zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)   ! (mol/L) * (m/s) 
     186            zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj) * zkgco2(ji,jj)   ! (mol/L) * (m/s) 
    185187            zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ? 
    186188            oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 
     
    189191 
    190192            ! Compute O2 flux  
    191             zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
     193            zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s) 
    192194            zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj) 
    193195            zoflx(ji,jj) = zfld16 - zflu16 
     
    222224         ENDIF 
    223225         IF( iom_use( "Dpco2" ) ) THEN 
    224            zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 
     226           zw2d(:,:) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1) 
    225227           CALL iom_put( "Dpco2" ,  zw2d ) 
    226228         ENDIF 
    227229         IF( iom_use( "Dpo2" ) )  THEN 
    228            zw2d(:,:) = ( atcox * patm(:,:) - trb(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1) 
     230           zw2d(:,:) = ( atcox * patm(:,:) - atcox * trn(:,:,1,jpoxy) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) 
    229231           CALL iom_put( "Dpo2"  , zw2d ) 
    230232         ENDIF 
     
    238240            trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1)  
    239241            trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1)  
    240             trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)  
     242            trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:) + rtrn ) ) * tmask(:,:,1)  
    241243         ENDIF 
    242244      ENDIF 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

    r6237 r6333  
    9191                  zalka = trb(ji,jj,jk,jptal) / zfact 
    9292                  ! CALCULATE [ALK]([CO3--], [HCO3-]) 
    93                   zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 
     93                  zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph / ( aphscale(ji,jj,jk) + rtrn )  & 
     94                  &       + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) ) 
    9495                  ! CALCULATE [H+] and [CO3--] 
    9596                  zaldi = zdic - zalk 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90

    r6237 r6333  
    7272      CHARACTER (len=25) :: charout 
    7373      REAL(wp), POINTER, DIMENSION(:,:  ) :: zpdep, zsidep, zwork1, zwork2, zwork3 
     74      REAL(wp), POINTER, DIMENSION(:,:)   :: zsedcal, zsedsi, zsedc 
    7475      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdenit2d, zironice, zbureff 
    7576      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsbio3, zwsbio4, zwscal 
     
    8384      ! Allocate temporary workspace 
    8485      CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
     86      CALL wrk_alloc( jpi, jpj, zsedcal,  zsedsi, zsedc ) 
    8587      CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
    8688      CALL wrk_alloc( jpi, jpj, jpk, zsoufer ) 
     
    9193      zwork2  (:,:) = 0.e0 
    9294      zwork3  (:,:) = 0.e0 
     95      zsedsi   (:,:) = 0.e0 
     96      zsedcal  (:,:) = 0.e0 
     97      zsedc    (:,:) = 0.e0 
    9398 
    9499      ! Iron input/uptake due to sea ice : Crude parameterization based on Lancelot et al. 
     
    298303            tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
    299304            tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
     305            zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss / zdep 
     306            zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss / zdep 
    300307#endif 
    301308         END DO 
     
    336343            tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) ) 
    337344            tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt 
    338             sdenit(ji,jj) = rdenit * zpdenit * fse3t(ji,jj,ikt) 
     345            sdenit(ji,jj) = rdenit * zpdenit / zdep 
     346            zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc / zdep 
    339347#endif 
    340348         END DO 
     
    392400               CALL iom_put( "INTNFIX" , zwork1 )  
    393401            ENDIF 
     402            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * 1.e+3 ) 
     403            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * 1.e+3 ) 
     404            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * 1.e+3 ) 
     405            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * 1.e+3 * rno3 ) 
    394406         ENDIF 
    395407      ELSE 
     
    405417      ! 
    406418      CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff ) 
     419      CALL wrk_dealloc( jpi, jpj, zsedcal , zsedsi, zsedc ) 
    407420      CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal ) 
    408421      CALL wrk_dealloc( jpi, jpj, jpk, zsoufer ) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90

    r5783 r6333  
    3838 
    3939   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget, po4budget 
    40    REAL(wp) :: xfact1, xfact2 
     40   REAL(wp) :: xfact1, xfact2, xfact3 
    4141   INTEGER ::  numco2, numnut, numnit  !: logical unit for co2 budget 
    4242 
     
    474474      !!--------------------------------------------------------------------- 
    475475      ! 
    476       INTEGER , INTENT( in ) ::   kt      ! ocean time-step index       
    477       REAL(wp)               ::  zfact        
    478       REAL(wp) ::  zrdenittot, zsdenittot, znitrpottot 
     476      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index       
     477      REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot 
    479478      CHARACTER(LEN=100)   ::   cltxt 
    480479      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol 
     
    492491            xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr 
    493492            xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr 
     493            xfact3 = 1.e+3 * rfact2r * rno3   ! conversion molC/l/kt ----> molN/m3/s 
    494494            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron' 
    495495            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt) 
     
    574574      IF( iom_use( "Sdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN 
    575575         zsdenittot   = glob_sum ( sdenit(:,:) * e1e2t(:,:) ) 
    576          CALL iom_put( "Sdenit", sdenit(:,:) * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments 
     576         CALL iom_put( "Sdenit", sdenit(:,:) * xfact3 * tmask(:,:,1) )  ! Nitrate reduction in the sediments 
    577577      ENDIF 
    578578 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90

    r5783 r6333  
    101101   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi         !: ??? 
    102102   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   excess     !: ??? 
     103   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aphscale   !:  
     104 
    103105 
    104106   !!* Temperature dependancy of SMS terms 
     
    154156         &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,       & 
    155157         &      akw3(jpi,jpj,jpk)    , borat (jpi,jpj,jpk) ,       & 
    156          &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,     STAT=ierr(4) ) 
     158         &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,       & 
     159         &      aphscale(jpi,jpj,jpk),                           STAT=ierr(4) ) 
    157160         ! 
    158161      !* Temperature dependancy of SMS terms 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90

    r5783 r6333  
    115115      po4r    =   1._wp / 122._wp 
    116116      o2nit   =  32._wp / 122._wp 
    117       rdenit  = 105._wp /  16._wp 
     117      o2ut    = 133._wp / 122._wp 
     118      rdenit  =  ( ( o2ut + o2nit ) * 0.80 - rno3 - rno3 * 0.60 ) / rno3 
    118119      rdenita =   3._wp /  5._wp 
    119       o2ut    = 133._wp / 122._wp 
     120 
    120121 
    121122      ! Initialization of tracer concentration in case of  no restart  
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90

    r5783 r6333  
    107107                
    108108               jl = n_trc_index(jn)  
    109                CALL trc_dta( kt, sf_trcdta(jl),rf_trfac(jl) )   ! read tracer data at nit000 
    110                ztrcdta(:,:,:) = sf_trcdta(jl)%fnow(:,:,:) 
     109               CALL trc_dta( kt, sf_trcdta(jl) )   ! read tracer data at nit000 
     110               ztrcdta(:,:,:) = sf_trcdta(jl)%fnow(:,:,:) * tmask(:,:,:) * rf_trfac(jl) 
    111111 
    112112               SELECT CASE ( nn_zdmp_tr ) 
     
    187187      INTEGER :: ji , jj, jk, jn, jl, jc                     ! dummy loop indicesa 
    188188      INTEGER :: isrow                                      ! local index 
    189       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta       ! 3D  workspace 
    190189 
    191190      !!---------------------------------------------------------------------- 
     
    278277         IF(lwp)  WRITE(numout,*) 
    279278         ! 
    280          CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )   ! Memory allocation 
    281          ! 
    282279         DO jn = 1, jptra 
    283280            IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file 
    284281                jl = n_trc_index(jn) 
    285                 CALL trc_dta( kt, sf_trcdta(jl),rf_trfac(jl) )   ! read tracer data at nit000 
    286                 ztrcdta(:,:,:) = sf_trcdta(jl)%fnow(:,:,:) 
     282                CALL trc_dta( kt, sf_trcdta(jl) )   ! read tracer data at nit000 
    287283                DO jc = 1, npncts 
    288284                   DO jk = 1, jpkm1 
    289285                      DO jj = nctsj1(jc), nctsj2(jc) 
    290286                         DO ji = nctsi1(jc), nctsi2(jc) 
    291                             trn(ji,jj,jk,jn) = ztrcdta(ji,jj,jk) * tmask(ji,jj,jk) 
     287                            trn(ji,jj,jk,jn) = sf_trcdta(jl)%fnow(ji,jj,jk) * tmask(ji,jj,jk) * rf_trfac(jl) 
    292288                            trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn) 
    293289                         ENDDO 
     
    297293             ENDIF 
    298294          ENDDO 
    299           CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta ) 
     295          ! 
    300296      ENDIF 
    301297      ! 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90

    r5781 r6333  
    5656      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    5757      !! 
    58       INTEGER            :: jn 
     58      INTEGER            :: ji, jj, jk, jn 
     59      REAL(wp)           :: zdep 
    5960      CHARACTER (len=22) :: charout 
    6061      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd 
     
    6667 
    6768      rldf = rldf_rat 
    68  
     69      ! 
     70      r_fact_lap(:,:,:) = 1. 
     71      DO jk= 1, jpk 
     72         DO jj = 1, jpj 
     73            DO ji = 1, jpi 
     74               IF( fsdept(ji,jj,jk) > 200. .AND. gphit(ji,jj) < 5. .AND. gphit(ji,jj) > -5. ) THEN 
     75                  zdep = MAX( fsdept(ji,jj,jk) - 1000., 0. ) / 1000. 
     76                  r_fact_lap(ji,jj,jk) = MAX( 1., rn_fact_lap * EXP( -zdep ) ) 
     77               ENDIF 
     78            END DO 
     79         END DO 
     80      END DO 
     81      ! 
    6982      IF( l_trdtrc )  THEN 
    7083         CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd ) 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/TRP/trcnam_trp.F90

    r5783 r6333  
    4040   REAL(wp), PUBLIC ::   rn_ahtrc_0          !: diffusivity coefficient for passive tracer (m2/s) 
    4141   REAL(wp), PUBLIC ::   rn_ahtrb_0          !: background diffusivity coefficient for passive tracer (m2/s) 
     42   REAL(wp), PUBLIC ::   rn_fact_lap         !: Enhanced zonal diffusivity coefficent in the equatorial domain 
    4243 
    4344   !                                        !!: ** Treatment of Negative concentrations ( nam_trcrad ) 
     
    7475      NAMELIST/namtrc_ldf/ ln_trcldf_lap  ,     & 
    7576         &                 ln_trcldf_bilap, ln_trcldf_level,     & 
    76          &                 ln_trcldf_hor  , ln_trcldf_iso  , rn_ahtrc_0, rn_ahtrb_0 
     77         &                 ln_trcldf_hor  , ln_trcldf_iso  , rn_ahtrc_0, rn_ahtrb_0,   & 
     78         &                 rn_fact_lap 
     79 
    7780      NAMELIST/namtrc_zdf/ ln_trczdf_exp  , nn_trczdf_exp 
    7881      NAMELIST/namtrc_rad/ ln_trcrad 
     
    127130         WRITE(numout,*) '      diffusivity coefficient                                 rn_ahtrc_0 = ', rn_ahtrc_0 
    128131         WRITE(numout,*) '      background hor. diffusivity                             rn_ahtrb_0 = ', rn_ahtrb_0 
     132         WRITE(numout,*) '      enhanced zonal diffusivity                             rn_fact_lap = ', rn_fact_lap 
    129133      ENDIF 
    130134 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90

    r5783 r6333  
    170170            END DO 
    171171         ENDIF 
     172         ! 
     173         CALL lbc_lnk( sbc_trc(:,:,jn), 'T', 1. ) 
    172174         !                                       Concentration dilution effect on tracers due to evaporation & precipitation  
    173175         DO jj = 2, jpj 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90

    r5783 r6333  
    6767         IF( lk_trabbl )        CALL trc_bbl( kstp )            ! advective (and/or diffusive) bottom boundary layer scheme 
    6868         IF( ln_trcdmp )        CALL trc_dmp( kstp )            ! internal damping trends 
    69          IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kstp )        ! internal damping trends on closed seas only 
    7069                                CALL trc_adv( kstp )            ! horizontal & vertical advection  
    7170                                CALL trc_ldf( kstp )            ! lateral mixing 
     
    7877                                CALL trc_nxt( kstp )            ! tracer fields at next time step      
    7978         IF( ln_trcrad )        CALL trc_rad( kstp )            ! Correct artificial negative concentrations 
     79         IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kstp )        ! internal damping trends on closed seas only 
    8080 
    8181#if defined key_agrif 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90

    r5781 r6333  
    116116   USE ldftra_oce , ONLY :  aeiw     =>   aeiw        !: eddy induced velocity coef. at w-points (m2/s)  
    117117   USE ldftra_oce , ONLY :  lk_traldf_eiv  =>  lk_traldf_eiv     !: eddy induced velocity flag 
     118   USE ldftra_oce , ONLY :  r_fact_lap     =>  r_fact_lap        !: enhanced zonal diffusivity coefficient 
    118119 
    119120   !* vertical diffusion * 
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/trcdta.F90

    r5783 r6333  
    151151 
    152152 
    153    SUBROUTINE trc_dta( kt, sf_dta, zrf_trfac ) 
     153   SUBROUTINE trc_dta( kt, sf_dta ) 
    154154      !!---------------------------------------------------------------------- 
    155155      !!                   ***  ROUTINE trc_dta  *** 
     
    165165      INTEGER                     , INTENT(in   ) ::   kt     ! ocean time-step 
    166166      TYPE(FLD), DIMENSION(1)   , INTENT(inout) ::   sf_dta     ! array of information on the field to read 
    167       REAL(wp)                  , INTENT(in   ) ::   zrf_trfac  ! multiplication factor 
    168167      ! 
    169168      INTEGER ::   ji, jj, jk, jl, jkk, ik    ! dummy loop indices 
     
    234233         ENDIF 
    235234         ! 
    236          sf_dta(1)%fnow(:,:,:) = sf_dta(1)%fnow(:,:,:) * zrf_trfac   !  multiplicative factor 
    237          ! 
    238235         IF( lwp .AND. kt == nit000 ) THEN 
    239236               clndta = TRIM( sf_dta(1)%clvar )  
  • branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r5783 r6333  
    6161      INTEGER ::   jk, jn, jl    ! dummy loop indices 
    6262      CHARACTER (len=25) :: charout 
    63       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrcdta   ! 4D  workspace 
    6463      !!--------------------------------------------------------------------- 
    6564      ! 
     
    121120        IF( ln_trcdta .AND. nb_trcdta > 0 ) THEN  ! Initialisation of tracer from a file that may also be used for damping 
    122121            ! 
    123             CALL wrk_alloc( jpi, jpj, jpk, ztrcdta )    ! Memory allocation 
    124             ! 
    125122            DO jn = 1, jptra 
    126123               IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file 
    127124                  jl = n_trc_index(jn)  
    128                   CALL trc_dta( nit000, sf_trcdta(jl),rf_trfac(jl) )   ! read tracer data at nit000 
    129                   ztrcdta(:,:,:) = sf_trcdta(jl)%fnow(:,:,:) 
    130                   trn(:,:,:,jn) = ztrcdta(:,:,:) * tmask(:,:,:)   
     125                  CALL trc_dta( nit000, sf_trcdta(jl) )   ! read tracer data at nit000 
     126                  trn(:,:,:,jn) = sf_trcdta(jl)%fnow(:,:,:) * tmask(:,:,:) * rf_trfac(jl) 
     127                  ! 
    131128                  IF( .NOT.ln_trcdmp .AND. .NOT.ln_trcdmp_clo ) THEN      !== deallocate data structure   ==! 
    132129                     !                                                    (data used only for initialisation) 
     
    138135               ENDIF 
    139136            ENDDO 
    140             CALL wrk_dealloc( jpi, jpj, jpk, ztrcdta ) 
     137            ! 
    141138        ENDIF 
    142139        ! 
Note: See TracChangeset for help on using the changeset viewer.