[MITgcm-support] heavy diffusion and water source
Patrick Rosendahl
Patrick.Rosendahl at zmaw.de
Sun Nov 2 14:40:16 EST 2008
Martin,
I have OLx=OLy=3. I monitored in integr_continuity() the variable
hDivFlow (which is based on vFld(j=1 and 2)) and it IS different
with/without hack.
The rise is flow dependent, at first slowly (the forcing hasnt reached
the entire domain), and then in the flow's steady state the rise is
constant.
I just compiled against the MITgcm_20081102 version (daily snapshot)
with the same result: the sealevel is rising.
I am going to live with this for now.
Patrick
Martin Losch wrote:
> 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
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
More information about the MITgcm-support
mailing list