[MITgcm-devel] woes with exf interpolating
Jean-Michel Campin
jmc at ocean.mit.edu
Thu Apr 8 10:03:42 EDT 2010
Hi Matt,
On Mon, Apr 05, 2010 at 04:22:12PM -0700, Matthew Mazloff wrote:
> Hi JMC,
>
> Thanks for your help! I had some trouble checking in SOSE but will try
> to revisit later this week.
>
> As for your questions, I do compile with package debug. I do not
> TARGET_NEC_SX , could that help? And yeah, xgOrigin = 0.0,
I was a little bit too quick: need also to put debugLevel=2 (in "data")
to check for valid range of lon & lat.
>
> I have not checked regarding optimization -- can I just NOOPTFILE
> exf_interp.F or are there other routines I should include?
Yes, you could try with just this one, since we suspect the problem
to come from there. But I tried a simple test (1-D in x, with exactly
the same _lon0 & _lon_inc for the input, and a regular dx=1/6 MITgcm
grid) and the cubic method gives me the right things. So, might need
to look more carefully to your set-up.
Also, can we be sure that the problem is already there before adding
optimisation-adjustment-correction (whatever you call it), i.e.,
at iteration zero of the optimisation ?
Cheers,
Jean-Michel
>
> And I can try your fix to see if anything changes...
>
> Thanks!
> -Matt
>
>
>
>
>
>
>
>
> On Apr 5, 2010, at 11:23 AM, Jean-Michel Campin wrote:
>
>> Hi Matt,
>>
>> I did take a look at exf_interp.F, and did not see any obvious
>> reason for a problem.
>> One thing that would help is to have your set-up checked-in,
>> (I remember Patrick creating this MITgcm_contrib/SOSE 2 months ago,
>> but it's still almost empty).
>> but in the mean time, could be useful to know the status of few
>> CPP options (TARGET_NEC_SX, ALLOW_DEBUG ...), and I assume
>> you are using the default xgOrigin=0 (in "data") ?
>>
>> Compiling the pkg "debug" would turn on '#define ALLOW_DEBUG'
>> and would check for valid range of long & lat in S/R exf_interp,
>> which could be useful.
>>
>> Also, did you check that with no compiler optimization you get
>> the same thing and same problem ? (there are few loop sequence,
>> with j loop inside an outer i loop that might be changed by the
>> compiler).
>>
>> There is also something you could try in exf_interp.F,
>> to put xG in interval [lon0-long_inc*0.5 , lon0-long_inc*0.5 + 360 [
>> instead of in [ lon_0 , lon_0+360 [
>> by changing lines 113-115, from:
>> 113 xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
>> 114 & + threeSixtyRS*2.
>> 115 xG(i,j,bi,bj) = lon_0+mod(xG(i,j,bi,bj),threeSixtyRS)
>> to the modified version:
>> 113 xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0+long_inc*0.5 _d 0
>> 114 & + threeSixtyRS*2.
>> 115 xG(i,j,bi,bj) = lon_0-long_inc*0.5 _d 0
>> 116 & + mod(xG(i,j,bi,bj),threeSixtyRS)
>> It will be more "symetric" (although I don't think it should be
>> a problem, but I could be wrong).
>>
>> Cheers,
>> Jean-Michel
>>
>> On Sat, Apr 03, 2010 at 12:22:47PM -0700, Matthew Mazloff wrote:
>>> Thanks Dimitris! I'll look for other explanations
>>>
>>> appreciate it,
>>> -Matt
>>>
>>>
>>> On Apr 3, 2010, at 12:14 PM, Menemenlis, Dimitris (3248) wrote:
>>>
>>>> Matt, exf_interp assumes periodicity in zonal direction of input
>>>> array
>>>> and provides 2 overlaps on either side, which should be sufficient
>>>> for
>>>> bicubic, so I don't understand where discontinuity comes from.
>>>>
>>>> C enlarge boundary
>>>> do j=1,ny_in
>>>> arrayin(0,j) = arrayin(nx_in,j)
>>>> arrayin(-1,j) = arrayin(nx_in-1,j)
>>>> arrayin(nx_in+1,j) = arrayin(1,j)
>>>> arrayin(nx_in+2,j) = arrayin(2,j)
>>>> enddo
>>>> do i=-1,nx_in+2
>>>> arrayin(i,0) = arrayin(i,1)
>>>> arrayin(i,-1) = arrayin(i,1)
>>>> arrayin(i,ny_in+1) = arrayin(i,ny_in)
>>>> arrayin(i,ny_in+2) = arrayin(i,ny_in)
>>>> enddo
>>>>
>>>> A simple test would be to add one more overlap on each side of
>>>> arrayin, and also to the related x_in and y_in, and see if it fixes
>>>> problem.
>>>>
>>>> Dimitris
>>>>
>>>> On Apr 3, 2010, at 10:45 AM, Matthew Mazloff wrote:
>>>>
>>>>> Hi JMC,
>>>>>
>>>>> I copied it below -- I think its pretty standard. Is there
>>>>> something
>>>>> I'm missing in prescribing the forcing dimensions -- do I need to
>>>>> manually provide an overlap, or provide overlap information?
>>>>>
>>>>> Thanks
>>>>> Matt
>>>>>
>>>>>
>>>>>
>>>>> login4% more data.exf
>>>>> # $Header: /u/gcmpack/MITgcm/verification/lab_sea/input/data.exf,v
>>>>> 1.7
>>>>> 2004/07/11 06:21:44 dimitri Exp
>>>>> $
>>>>> #
>>>>> # *********************
>>>>> # External Forcing Data
>>>>> # *********************
>>>>> &EXF_NML_01
>>>>> #
>>>>> useExfYearlyFields = .TRUE.,
>>>>> useExfCheckRange = .FALSE.,
>>>>> exf_iprec = 32,
>>>>> exf_yftype = 'RL',
>>>>> useRelativeWind = .TRUE.,
>>>>> /
>>>>> #
>>>>> &EXF_NML_02
>>>>> #
>>>>> atempstartdate1 = 20050101,
>>>>> atempstartdate2 = 00000,
>>>>> atempperiod = 21600.0,
>>>>> #
>>>>> aqhstartdate1 = 20050101,
>>>>> aqhstartdate2 = 00000,
>>>>> aqhperiod = 21600.0,
>>>>> #
>>>>> precipstartdate1 = 20050101,
>>>>> precipstartdate2 = 00000,
>>>>> precipperiod = 21600.0,
>>>>> #
>>>>> uwindstartdate1 = 20050101,
>>>>> uwindstartdate2 = 00000,
>>>>> uwindperiod = 21600.0,
>>>>> #
>>>>> vwindstartdate1 = 20050101,
>>>>> vwindstartdate2 = 00000,
>>>>> vwindperiod = 21600.0,
>>>>> #
>>>>> runoffstartdate1 = 20050101,
>>>>> runoffstartdate2 = 00000,
>>>>> runoffperiod = 0.0,
>>>>> #
>>>>> swdownstartdate1 = 20050101,
>>>>> swdownstartdate2 = 00000,
>>>>> swdownperiod = 21600.0,
>>>>> #
>>>>> lwdownstartdate1 = 20050101,
>>>>> lwdownstartdate2 = 00000,
>>>>> lwdownperiod = 21600.0,
>>>>> #
>>>>> atempfile = 'SO6input/NCEP_R6_tmp2m_degC_SO',
>>>>> aqhfile = 'SO6input/NCEP_R6_spfh2m_SO',
>>>>> uwindfile = 'SO6input/NCEP_R6_u10m_SO',
>>>>> vwindfile = 'SO6input/NCEP_R6_v10m_SO',
>>>>> precipfile = 'SO6input/NCEP_R6_rain_SO_DT',
>>>>> runoffFile = 'SO6input/run-off_SO.6',
>>>>> lwdownfile = 'SO6input/NCEP_R6_dlw_SO',
>>>>> swdownfile = 'SO6input/NCEP_R6_dsw_SO',
>>>>> #
>>>>> /
>>>>> #
>>>>> &EXF_NML_03
>>>>> #
>>>>> exf_inscal_uwind = 1.,
>>>>> exf_inscal_vwind = 1.,
>>>>> exf_inscal_runoff = 3.170979198E-8,
>>>>> exf_offset_atemp = 273.15,
>>>>> #
>>>>> /
>>>>> #
>>>>> &EXF_NML_04
>>>>> #
>>>>> atemp_lon0 = 0.5D0,
>>>>> atemp_lon_inc = 1.D0,
>>>>> atemp_lat0 = -79.5D0,
>>>>> atemp_lat_inc = 60*1.D0,
>>>>> atemp_nlon = 360,
>>>>> atemp_nlat = 60,
>>>>> #
>>>>> aqh_lon0 = 0.5D0,
>>>>> aqh_lon_inc = 1.D0,
>>>>> aqh_lat0 = -79.5D0,
>>>>> aqh_lat_inc = 60*1.D0,
>>>>> aqh_nlon = 360,
>>>>> aqh_nlat = 60,
>>>>> #
>>>>> precip_lon0 = 0.5D0,
>>>>> precip_lon_inc = 1.D0,
>>>>> precip_lat0 = -79.5D0,
>>>>> precip_lat_inc = 60*1.D0,
>>>>> precip_nlon = 360,
>>>>> precip_nlat = 60,
>>>>> #
>>>>> uwind_lon0 = 0.5D0,
>>>>> uwind_lon_inc = 1.D0,
>>>>> uwind_lat0 = -79.5D0,
>>>>> uwind_lat_inc = 60*1.D0,
>>>>> uwind_nlon = 360,
>>>>> uwind_nlat = 60,
>>>>> #
>>>>> vwind_lon0 = 0.5D0,
>>>>> vwind_lon_inc = 1.D0,
>>>>> vwind_lat0 = -79.5D0,
>>>>> vwind_lat_inc = 60*1.D0,
>>>>> vwind_nlon = 360,
>>>>> vwind_nlat = 60,
>>>>> #
>>>>> swdown_lon0 = 0.5D0,
>>>>> swdown_lon_inc = 1.D0,
>>>>> swdown_lat0 = -79.5D0,
>>>>> swdown_lat_inc = 60*1.D0,
>>>>> swdown_nlon = 360,
>>>>> swdown_nlat = 60,
>>>>> #
>>>>> lwdown_lon0 = 0.5D0,
>>>>> lwdown_lon_inc = 1.D0,
>>>>> lwdown_lat0 = -79.5D0,
>>>>> lwdown_lat_inc = 60*1.D0,
>>>>> lwdown_nlon = 360,
>>>>> lwdown_nlat = 60,
>>>>> #
>>>>> /
>>>>> #
>>>>> &EXF_NML_OBCS
>>>>> #
>>>>> obcsNstartdate1 = 20050101,
>>>>> obcsNstartdate2 = 00000,
>>>>> obcsNperiod = 2635200.0,
>>>>> #
>>>>> /
>>>>>
>>>>>
>>>>>
>>>>> On Apr 3, 2010, at 7:07 AM, Jean-Michel Campin wrote:
>>>>>
>>>>>> Hi Matt,
>>>>>>
>>>>>> Could you send also the "data.exf" you are using ? or could I find
>>>>>> it
>>>>>> somewhere ?
>>>>>> Thanks,
>>>>>> Jean-Michel
>>>>>>
>>>>>> On Fri, Apr 02, 2010 at 03:17:51PM -0700, Matthew Mazloff wrote:
>>>>>>> Hello,
>>>>>>>
>>>>>>> A problem has been found (thanks Fabien!) regarding a
>>>>>>> discontinuity
>>>>>>> in
>>>>>>> the wind fields around longitude 0/360.
>>>>>>>
>>>>>>> My set-up is 1/6 of a degree and I use NCEP on a 1 degree
>>>>>>> grid and
>>>>>>> use
>>>>>>> exf interpolation. Within ~1 degree east and west of longitude 0
>>>>>>> there
>>>>>>> is a discontinuity in the wind field -- it is small but quite
>>>>>>> significant in the x-derivative. My guess is that the
>>>>>>> interpolation
>>>>>>> routine is not properly accounting for periodicity -- and instead
>>>>>>> perhaps incorporating zeros.
>>>>>>>
>>>>>>> Some important stuff -- I am using NCEP winds -- not stress --
>>>>>>> and I
>>>>>>> have switched to interp_method=2 instead of the hard-coded
>>>>>>> interp_method=1. This is because with interp_method 1 (linear)
>>>>>>> there
>>>>>>> are discontinuities in the wind-stress curl everywhere. With
>>>>>>> bicubic
>>>>>>> interp the curl is smooth -- with this exception around
>>>>>>> longitude 0
>>>>>>>
>>>>>>> I was hoping someone with more expertise could shed light on this
>>>>>>> before
>>>>>>> I go digging in the code
>>>>>>>
>>>>>>> Thanks!
>>>>>>> -Matt
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> MITgcm-devel mailing list
>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>
>>>>>> _______________________________________________
>>>>>> MITgcm-devel mailing list
>>>>>> MITgcm-devel at mitgcm.org
>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> MITgcm-devel mailing list
>>>>> MITgcm-devel at mitgcm.org
>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>
>>>>
>>>> _______________________________________________
>>>> MITgcm-devel mailing list
>>>> MITgcm-devel at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>
>>>
>>> _______________________________________________
>>> MITgcm-devel mailing list
>>> MITgcm-devel at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list