[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