[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