[MITgcm-support] Re: Re.: hFac's and obcs

[Thomas Haine]

Hi Martin, 

Turning off rStar solves the NaNs in the initial W dump. So, they must
be caused in the rStar code. However, the cg2d solver then diverges at
its first call and W is NaN throughout the domain. 

So, changing the hFacS this way causes severe issues in at least 2
separate parts of the code. 

My expertise in these topics is feeble so perhaps someone else can
comment?  Alistair?

Tom.

On Wed, 2005-09-14 at 09:09 +0200, Martin Losch wrote:
> Hi Tom,
> 
> I have no further ideas, but between initialization and the first 
> monitor, w is computed from the initial horizontal velocity field:
> in the_main_loop:
> call initialise_varia
> call monitor
> call do_the_model_io
> 
> in initialise_varia:
> ini_dynvars
> packages_init_variables -> obcs_init_variables
> (convective_adjustment, if you use that)
> r*-business
> integr_continuity -> calculation of w
> 
> For debugging, I personally would turn off the r*-stuff, to make things 
> easier (and because i still don't understand what's going on there 
> exactly, but you probably do). The diverging w's can mean that your 
> u,v-fields have strong divergences on the boundary, or that the 
> topography has a gradient normal to the boundary. Or the diverging w's 
> can certainly mean, that my idea about the hfacs's was complete 
> nonsense.
> 
> Martin
> 
> On Sep 13, 2005, at 10:04 PM, Thomas Haine wrote:
> 
> > Hi Martin,
> >
> > Thanks for your comments. Modifying hFacS in the way you suggest causes
> > the vertical velocity field to diverge over most of the southern
> > boundary at step 0 (i.e. it's visible in the MONITOR output at step 0
> > and in the 1st dump).
> >
> > I tried this variation on your theme:
> >
> > hFacS(I,1,K,bi,bj)=hFacC(I,1,K,bi,bj)
> >
> > but got the same results. Not much goes on between initialising and
> > calculating the MONITOR but I've not yet found where W gets set wrong 
> > or
> > understood why.
> >
> > Any other ideas?
> >
> > Thanks, Tom.
> >
> >
> >> -------- Forwarded Message --------
> >>> From: Martin Losch <mlosch at awi-bremerhaven.de>
> >>> Reply-To: mitgcm-support at mitgcm.org
> >>> To: mitgcm-support at mitgcm.org
> >>> Subject: Re: [MITgcm-support] hFac's and obcs
> >>> Date: Mon, 12 Sep 2005 10:45:55 +0200
> >>>
> >>> Hi Tom,
> >>>
> >>> I tried to find places in the OBCS-related code, where hFacS (or 
> >>> maskS)
> >>> could be a problem. Having had only a quick look I couldn't find
> >>> anything. In fact, the open southern boundary values for v-velocity 
> >>> are
> >>> prescribed at OB_Js+1, that is at J=2 in your case, so that the value
> >>> of hFacS/maskS at J=1 shouldn't matter:
> >>> In obcs_apply_uv.F, the velocities fields are set:
> >>>>        IF (OB_Js(I,bi,bj).NE.0) THEN
> >>>>         vFld(I,OB_Js(I,bi,bj)+1,K,bi,bj)=OBSv(I,K,bi,bj)
> >>>>      &                              
> >>>> *_maskS(I,OB_Js(I,bi,bj)+1,K,bi,bj)
> >>>>         uFld(I,OB_Js(I,bi,bj),K,bi,bj)=OBSu(I,K,bi,bj)
> >>>>      &                              
> >>>> *_maskW(I,OB_Js(I,bi,bj),K,bi,bj)
> >>>>         vFld(I,OB_Js(I,bi,bj),K,bi,bj)=OBSv(I,K,bi,bj)
> >>>>      &                              
> >>>> *_maskS(I,OB_Js(I,bi,bj)+1,K,bi,bj)
> >>>>        ENDIF
> >>> In your case, OB_Js = 1 for all i, correct?
> >
> > Yes, that's right.
> >
> >>> That means, that
> >>> vFld(:,1,:)=vFld(:,2,:)=OBSv(:,:).*maskS(:,2,:), so that 
> >>> maskS(:,1,:),
> >>> which should really be 0 for most i, is never used.]
> >>> In solve_for_pressure.F, cg2d_b and cg2d_x are set to zero on the
> >>> boundaries and hFacS/W are not used.
> >>>
> >>> However, both in mom_fluxform and mom_vecinv, hfacs/hfacw are
> >>> excessively used, so that I would agree, there is a problem with 
> >>> using
> >>> OBCS. Have you tried to modify the hFac's in ini_masks_etc? After the
> >>> line
> >>> CALL EXCH_UV_XYZ_RS(hFacC,hFacS,.FALSE.,myThid)
> >>> you could hack your hfacs according to your domain. Something like:
> >>>>       DO bj=myByLo(myThid), myByHi(myThid)
> >>>>        DO bi=myBxLo(myThid), myBxHi(myThid)
> >>>>         DO K=1, Nr
> >>>>          DO J=1-Oly,1
> >>>>           DO I=1-Olx,sNx+Olx
> >>>>            hFacS(I,J,K,bi,bj)=hFacS(I,2,K,bi,bj)
> >>>>           ENDDO
> >>>>          ENDDO
> >>>>          DO J=1-Oly,sNy+Oly
> >>>>           DO I=1-Olx,1
> >>>>            hFacW(I,J,K,bi,bj)=hFacW(2,J,K,bi,bj)
> >>>>           ENDDO
> >>>>          ENDDO
> >>>>         ENDDO
> >>>>        ENDDO
> >>>>       ENDDO
> >>> If the code survives such a such a fundamental modification, that is,
> >>> if this fixes the problem, we should consider including it in the
> >>> OBCS-package.
> >>>
> >>> Martin
-- 
--------------------------------------------
 Thomas W. N. Haine
                                                                                                                                                              
 Associate Professor of Physical Oceanography,
 Department of Earth & Planetary Sciences,
 329 Olin Hall, The Johns Hopkins University,
 Baltimore, MD 21218, USA.
 Tel : 410 516 7048,   Fax : 410 516 7933
 Thomas.Haine at jhu.edu
 http://www.jhu.edu/~eps/faculty/haine
--------------------------------------------