[MITgcm-devel] woes with exf interpolating
Jean-Michel Campin
jmc at ocean.mit.edu
Mon Apr 5 14:23:09 EDT 2010
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
More information about the MITgcm-devel
mailing list