[MITgcm-support] heavy diffusion and water source
Martin Losch
Martin.Losch at awi.de
Sun Nov 2 13:25:36 EST 2008
Patrick,
I have no clue what is going on, you should not have to do this hack
about the v(2)=v(1), it should be like that anyway, please test that
(vfld(:,sNy+1:sNy+Oly,:,:,:) = vFld(:,1:Oly,:,:,:), i.e. your case
sNy=1 and OLy=3(?, otherwise advScheme 33 should not work), vfld(:,
2:4,:,:,:)=vfld(:,1:3,:,:,:) (at least after an exchange).
2.5e-5 in 40000 sec means 6.25e-9 in one time step (deltaT=10sec),
that seems to be too much, doesn't it? It should be on the order of
1e-14 for a double precision code. Is the rise linear in time, or
flow dependent?
Martin
On 2 Nov 2008, at 18:37, Patrick Rosendahl wrote:
> Hello Martin,
> thanks for the help.
>
> I am working with periodic boundaries.
>
> A) I have added useRealFreshWaterFlux, select_rStar, nonlinFreeSurf
> and hFac* to the data file (the rest was already set as suggested,
> also in my CPP_OPTIONS.h). But I still get +2.5e-5 m sealevel rise
> after t=40000s.
> I guess it is because I am using a 1-pixel model in the XZ plane.
> When I am forcing vVel(1) = vVel(2) everything is good:
>
> -------------------------------------integr_continuity.F
> #ifdef EXACT_CONSERV
> IF (exactConserv) THEN
>
> C-- Compute the Divergence of The Barotropic Flow :
>
> C- force duplicated values, because we are working with a slice
> only
> IF ( sNy.eq.1 ) THEN
> DO k=1,Nr
> DO i=1,sNx+1
> vFld(i,2,k,bi,bj) = vFld(i,1,k,bi,bj)
> ENDDO
> ENDDO
> ENDIF
> -------------------------------------integr_continuity.F
>
> However, still wondering why the level rises... Maybe the level is
> computed at another point in the code with a different formula.
> What I mean is this situation:
> "The formula" might be partially evaluated before and after the
> model duplicates the vVel-values itself. Then the values are
> inconsistent and the model "thinks" there was a rise before the
> timestep already - hence it conserves it.
>
>
> B) It turned out that the "diffusivity" was rather fluxes due to
> rising sea level.
>
> Cheers,
> Patrick
>
> Martin Losch wrote:
>> Patrick,
>> not sure what's happening, but according to your data.pkg you are
>> NOT using OBCS (so no "open boundaries"), probably you mean
>> periodic boudaries?
>> A) exact conservation of volume (no sea level rise for a closed
>> system) is achieved like this:
>> 1. in CPP_OPTIONS.h: #define EXACT_CONSERV (that's the default)
>> 2. in data: exactConserv = .true. (not default for backward
>> compatibility, I guess), and
>> useRealFreshWaterFlux=.TRUE., if you have surface fresh water
>> forcing (salinity restoring or EmPmR)
>> if that's not sufficient, then you could try the non-linear free
>> surface (may not really exactly work with OBCS, if that's required):
>> 3. in CPP_OPTIONS.h: #define NONLIN_FRSURF (not default)
>> 4. in data (taken from verification/global_ocean.cs32x15/input/
>> data: select_rStar=2, nonlinFreeSurf=4, hFacInf=0.2, hFacSup=2.0,
>> please look up the meaning of the parameters in the documentation
>> B) AdvScheme=2 is the simplest scheme (2nd order central
>> differences) and should not be stable without additional
>> diffusion, but that's not your problem, I guess. You use 33 for
>> temp and saltAdvScheme, this scheme should be used with
>> staggerTimeStep=.true., otherwise you'll get a lot of vertical
>> mixing (by static instability?), often discussed in the email-
>> list, but never quite understood (at least not by myself), this
>> mixing will affect your ptracers regardless of their advection
>> scheme. Setting staggerTimeStep=.true. should fix your problem.
>> I'd then go with 33 for ptracers as well; there is numerical
>> diffusion in that scheme, but it's small and directed towards
>> stabilizing the scheme.
>> Martin
>> On 2 Nov 2008, at 02:38, Patrick Rosendahl wrote:
>>> Hello,
>>>
>>> I have a flow in wind direction (constant tx) with open
>>> boundaries. And I have 2 problems:
>>>
>>> A) The sealevel rises very slowly like 1e-6 in 100s. I fixed this
>>> by setting hDivFlow(i,j)=0.0 and same in integrate_for_w(). This
>>> kept my eta=w=0. But this is not a clean method. V is zero as it
>>> should. Any thoughts on this?
>>>
>>> B) There is a lot of vertical diffusivity in the ptracers,
>>> although the transport is only in x-direction. I have set the
>>> diffusivity to 0.0 - no help. I am using advectionScheme=2, since
>>> it appears to have the least diffusivity. Why is there so much
>>> diffusivity?
>>>
>>> Thanks,
>>> Patrick
>>>
>>> --------------------------------------------------------------------
>>> ---------------------
>>> DATA:
>>>
>>> # ====================
>>> # | Model parameters |
>>> # ====================
>>> #
>>> # Continuous equation parameters
>>> &PARM01
>>> tRef=51*20.0,
>>> sRef=51*0.0,
>>> viscAh=1.3E-6,
>>> viscAz=1.3E-6,
>>> no_slip_sides=.FALSE.,
>>> no_slip_bottom=.TRUE.,
>>> diffKhT=1.4E-7,
>>> diffKzT=1.4E-7,
>>> diffKhS=1.1E-9,
>>> diffKzS=1.1E-9,
>>> f0=0.0,
>>> sBeta =0.,
>>> gravity=9.81,
>>> rhoConst=1000.0,
>>> rhoNil=1000.0,
>>> # heatCapacity_Cp=3900.0,
>>> implicitDiffusion=.TRUE.,
>>> implicitViscosity=.TRUE.,
>>> eosType='LINEAR',
>>> readBinaryPrec=32,
>>> # hFacMin=0.2,
>>> # hFacInf=0.2,
>>> tempAdvScheme=33,
>>> saltAdvScheme=33,
>>>
>>> nonHydrostatic=.FALSE.,
>>> # nonlinFreeSurf=3,
>>> exactConserv=.TRUE.,
>>> implicitFreeSurface=.TRUE.,
>>> rigidLid=.FALSE.,
>>> &
>>>
>>> # Elliptic solver parameters
>>> &PARM02
>>> cg2dMaxIters=200,
>>> cg2dTargetResidual=1.E-13,
>>> cg3dMaxIters=10,
>>> cg3dTargetResidual=1.E-9,
>>> &
>>>
>>> # Time stepping parameters
>>> &PARM03
>>> nIter0=0,
>>> nTimeSteps=1000,
>>> deltaT=10.0,
>>> abEps=0.1,
>>> pChkptFreq=500000.0,
>>> chkptFreq=500000.0,
>>> dumpFreq=100.0,
>>> monitorFreq=500000.0,
>>> # outputTypesInclusive=.TRUE.,
>>> &
>>>
>>> # Gridding parameters
>>> &PARM04
>>> usingCartesianGrid=.TRUE.,
>>> usingCylindricalGrid=.FALSE.,
>>> usingCurvilinearGrid=.FALSE.,
>>> delX=20*20.0,
>>> delY=1*20.0,
>>> delZ=51*0.2,
>>> &
>>>
>>> # Input datasets
>>> &PARM05
>>> zonalWindFile='windstressx.const_1.000.bin.real*4',
>>> bathyFile='depth.open_10.0.bin.real*4',
>>> &
>>>
>>> --------------------------------------------------------------------
>>> ---------------------
>>> DATA.PTRACERS:
>>>
>>> &PTRACERS_PARM01
>>> # PTRACERS_dumpFreq,
>>> # PTRACERS_taveFreq,
>>> # PTRACERS_monitorFreq,
>>> PTRACERS_advScheme(1)=2,
>>> PTRACERS_advScheme(2)=2,
>>> # PTRACERS_ImplVertAdv,
>>> PTRACERS_diffKh(1)=0.0,
>>> PTRACERS_diffKh(2)=0.0,
>>> PTRACERS_diffK4(1)=0.0,
>>> PTRACERS_diffK4(2)=0.0,
>>> PTRACERS_diffKr(1)=0.0,
>>> PTRACERS_diffKr(2)=0.0,
>>> # PTRACERS_ref,
>>> # PTRACERS_useGMRedi=.FALSE.,
>>> # PTRACERS_useKPP=.FALSE.,
>>> # PTRACERS_Iter0,
>>> PTRACERS_numInUse=2,
>>> # PTRACERS_initialFile,
>>> # PTRACERS_useRecords,
>>> PTRACERS_names(1)='tracer1',
>>> PTRACERS_names(2)='tracer2',
>>> PTRACERS_long_names(1)='tracer1',
>>> PTRACERS_long_names(2)='tracer2',
>>> PTRACERS_units(1)='ppm',
>>> PTRACERS_units(2)='ppm',
>>> # PTRACERS_timeave_mnc,
>>> PTRACERS_snapshot_mnc=.TRUE.,
>>> PTRACERS_monitor_mnc=.TRUE.,
>>> PTRACERS_pickup_write_mnc=.TRUE.,
>>> # & PTRACERS_pickup_read_mnc
>>> &
>>>
>>> --------------------------------------------------------------------
>>> ---------------------
>>> DATA.PKG:
>>>
>>> # Packages
>>> &PACKAGES
>>> useMNC=.TRUE.,
>>> # useOBCS=.TRUE.,
>>> # useKPP=.TRUE.,
>>> # useMY82=.TRUE.,
>>> usePTRACERS=.TRUE.,
>>> &
>>> _______________________________________________
>>> MITgcm-support mailing list
>>> MITgcm-support at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-support
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
More information about the MITgcm-support
mailing list