Работа с GrADS

Метеоцентр.Азия - наш новый сайт с высокодетализированными прогнозами погоды по пунктам России и мира
Облегчённая версия Метеоклуба (для КПК)
Карта активных участников Метеоклуба (105 кБ)    Таблица дней рождения активных участников Метеоклуба
Клуб любителей метеорологии (группа ВКонтакте)

Работа с GrADS

Сейчас в Метеоклубе:
Участников - 2 [ Adrenalines, Balin ]
Максимальное одновременное количество посетителей: 3 [20 Янв 2017 08:32]
Гостей - 0 / Участников - 3

 - Начало - Ответить - Статистика - Регистрация - Поиск -
МЕТЕОКЛУБ : независимое сообщество любителей метеорологии (Европа и Азия) : ФОРУМ О ПОГОДЕ И ПРИРОДЕ / Компьютерная техника и интернет в метеорологии / Работа с GrADS
<< . 1 . 2 . 3 . 4 . 5 . 6 . 7 . 8 . 9 . 10 . >>
Автор Сообщение
bullterrier
Участник
Письмо
Пермь
# Дата: 20 Ноя 2012 13:57


Скажите, а в OpenGrads есть фунцкция, которая бы по температуре и точке росы выдавала влажность в процентах? (Канадская модель GEM). Покопался в документации здесь: http://www.iges.org/grads/gadoc/gadocindex.html , но ничего не нашел, возможно из-за плохого знания английского.

Cumulonimbus incus
Участник
Письмо
Кишинёв, Молдова
# Дата: 20 Ноя 2012 17:16


bullterrier

Предлагаю такой способ. Здесь t и td - темп-ра и точка росы в градусах Цельсия.

'define rh = 100*exp((18.678-td/234.5)*(td/(257.14+td))-(18.678-t/234.5)*(t/(257.14+t)))'

bullterrier
Участник
Письмо
Пермь
# Дата: 21 Ноя 2012 15:50 - Поправил: bullterrier


Cumulonimbus incus

Спасибо, работает!

Раз этой функции у Grads нет, то можно сделать так (я сделал это на основе функции dewpt)

function humid (tk, tdk)

if (tdk='tdk'|tdk='' )
say 'Purpose: humidity [%]'
say 'Usage: display humid(t,rh)'
say ' tk = temperature [K]'
say ' tdk = dewpoint temperature [K]'
return
else
tdkt=tdk
endif

'define tc = 'tk' - 273.15'
'define tdc = 'tdkt' - 273.15'
'define rhp = 100*exp((18.678-tdc/234.5)*(tdc/(257.14+tdc))-(18. 678-tc/234.5)*(tc/(257.14+tc)))'


'undefine tc'
'undefine tdc'

return 'rhp'

Затем сохранить ее в файл humid.gsf и поместить в папку со скриптами (/opengrads/Resources/Scripts)

PS. Так не получится схитрить, к сожалению. OpenGrads вылетает

bullterrier
Участник
Письмо
Пермь
# Дата: 22 Ноя 2012 20:31


Подсчет Surface-Based LI в OpenGrads.



<...>

'define Temp500V='virtual2('t500mb+273.15','dp500mb+273.15 ',500)'-273.15'
'define TLcl='Templcl('tsfc','dpsfc')
'define PLcl='Preslcl('tsfc','dpsfc','psfc')

'define PclTemp='LiftWet('TLcl','PLcl',500)
'define PclTempV='virtual2('PclTemp+273.15','PclTemp+273.1 5',500)'-273.15'
'define LIMax=Temp500V-PclTempV'

'define variable=LIMax'

<...>

function LiftWet(startt,startp,endp)

*------------------------------------------------- -------------------
* Lift a parcel moist adiabatically from startp to endp.
* Init temp is startt in C. If you wish to see the parcel's
* path plotted, display should be 1. Returns temp of parcel at endp.
*------------------------------------------------- -------------------

temp=startt
pres=startp
delp=10
While (pres >= endp)
temp=temp-100*delp*gammaw(temp,pres-delp/2,100)
pres=pres-delp
EndWhile
return(temp)

function gammaw(tempc,pres,rh)

*------------------------------------------------- ----------------------
* Function to calculate the moist adiabatic lapse rate (deg C/Pa) based
* on the temperature, pressure, and rh of the environment.
*------------------------------------------------- ---------------------

tempk=tempc+273.15
es=satvap2(tempc)
ws=mixratio(es,pres)
w=rh*ws/100
tempv=virtual(tempk,w)
latent=latentc(tempc)

A=1.0+latent*ws/(287*tempk)
B=1.0+0.622*latent*latent*ws/(1005*287*tempk*tempk )
Density=100*pres/(287*tempv)
lapse=(A/B)/(1005*Density)
return(lapse)

function latentc(tempc)

*------------------------------------------------- ----------------------
* Function to return the latent heat of condensation in J/kg given
* temperature in degrees Celsius.
*------------------------------------------------- ----------------------

val=2502.2-2.43089*tempc*1000

return(val)

function Templcl(temp,dewp)

*------------------------------------------------- -----
* Calculate the temp at the LCL given temp & dewp in C
*------------------------------------------------- -----

tempk=temp+273.15
dewpk=dewp+273.15
Parta=1/(dewpk-56)
Partb=log(tempk/dewpk)/800
Tlcl=(1/(Parta+Partb)+56)-273.15
return(Tlcl)

function Preslcl(temp,dewp,pres)

*------------------------------------------------- ------
* Calculate press of LCL given temp & dewp in C and pressure
*------------------------------------------------- ------

Tlclk=Templcl(temp,dewp)+273.15
tempk=temp+273.15
theta=tempk*pow(1000/pres,0.286)
plcl=1000*pow(Tlclk/theta,3.48)
return(plcl)

function mixratio(e,p)

*------------------------------------------------- -----
* Given vapor pressure and pressure, return mixing ratio
*------------------------------------------------- ------

mix=0.622*e/(p-e)

return(mix)

function virtual2(temp,dewp,pres)

*------------------------------------------------- -----------
* Function to return virtual temperature in kelvin given temperature in
* kelvin and dewpoint in kelvin and pressure in mb
*------------------------------------------------- ------------

if (temp > 0 & dewp > 0)
* vap=satvap2(dewp-273.15)*
mix=mixratio(vap,pres)
tempv=virtual(temp,mix)
else
tempv=-9999
endif

return (tempv)

function virtual(temp,mix)

*------------------------------------------------- -----------
* Function to return virtual temperature given temperature in
* kelvin and mixing ratio in g/g.
*------------------------------------------------- ------------

tempv=temp*(1.0+0.6*mix)

return (tempv)

function mixratio(e,p)

*------------------------------------------------- -----
* Given vapor pressure and pressure, return mixing ratio
*------------------------------------------------- ------

mix=0.622*e/(p-e)

return(mix)

function satvap2(temp)

*------------------------------------------------- --------------
* Given temp in Celsius, returns saturation vapor pressure in mb
*---------------------------------------------------------------

es=6.112*exp(17.67*temp/(temp+243.5))

return(es)


http://gradsusr.org/pipermail/gradsusr/2012-March/ 015320.html

bullterrier
Участник
Письмо
Пермь
# Дата: 22 Ноя 2012 20:34


Попытался прикрутить этот код к скрипту для модели GEM, результыты получаются неверные.

<...>

'vtmp500 = 'virtual2('tmp500','dew500',500)'-273.15'

'tl=tlcl(t2m+273.15,rh2m)-273.15'
'pl=plcl(t2m+273.15,rh2m,spr)'

'PclTemp='LiftWet('tl','pl',500)
'PclTempV='virtual2('PclTemp+273.15','PclTemp+273. 15',500)'-273.15'
'sli=vtmp500-PclTempV'

<...>

function LiftWet(startt,startp,endp)

*------------------------------------------------- -------------------
* Lift a parcel moist adiabatically from startp to endp.
* Init temp is startt in C. If you wish to see the parcel's
* path plotted, display should be 1. Returns temp of parcel at endp.
*------------------------------------------------- -------------------

'temp='startt
'pres='startp
'delp=10'
'eendp='endp
i=0
While (i=0)
'gw='gammaw('temp','pres-delp/2',100)
'temp=temp-100*delp*gw'
'pres=pres-delp'
if ('pres<eendp')
i=1
endif
EndWhile
return('temp')

function Templcl(temp,dewp)

*------------------------------------------------- -----
* Calculate the temp at the LCL given temp & dewp in C
*------------------------------------------------- -----

'tempk='temp'+273.15'
'dewpk='dewp'+273.15'
'Parta=1/(dewpk-56)'
'Partb=log(tempk/dewpk)/800'
'Tlcl=1/(Parta+Partb)+56'
return('Tlcl-273.15')

function Preslcl(temp,dewp,pres)

*------------------------------------------------- ------
* Calculate press of LCL given temp & dewp in C and pressure
*------------------------------------------------- ------

'Tlclk='Templcl(temp,dewp)'+273.15'
'tempk='temp'+273.15'
'theta=tempk*pow(1000/'pres',0.286)'
'plcl=1000*pow(Tlclk/theta,3.48)'
return('plcl')

function gammaw(tempc,pres,rh)

*------------------------------------------------- ----------------------
* Function to calculate the moist adiabatic lapse rate (deg C/Pa) based
* on the temperature, pressure, and rh of the environment.
*------------------------------------------------- ---------------------

*'tempk='tempc'+273.15'
*'es='satvap2(tempc)
*'ws='mixratio('es',pres)
*'w='rh'*ws/100'
*'tempv='virtual('tempk','w')
*'latent='latentc(tempc)
*'A=1.0+latent*ws/(287*tempk)'
*'B=1.0+0.622*latent*latent*ws/(1005*287*tempk*tem pk)'
*'Density=100*'pres'/(287*tempv)'
*'lapse=(A/B)/(1005*Density)'

'ttempc='tempc
'ppres='pres
'rrh='rh
'tempk=ttempc+273.15'
'es='satvap2('ttempc')
'ws='mixratio('es','ppres')
'w=rrh*ws/100'
'tempv='virtual('tempk','w')
'latent='latentc('ttempc')
'A=1.0+latent*ws/(287*tempk)'
'B=1.0+0.622*latent*latent*ws/(1005*287*tempk*temp k)'
'Density=100*ppres/(287*tempv)'
'lapse=(A/B)/(1005*Density)'

return('lapse')

function latentc(tempc)

*------------------------------------------------- ----------------------
* Function to return the latent heat of condensation in J/kg given
* temperature in degrees Celsius.
*------------------------------------------------- ----------------------

*'val=2502.2-2.43089*'tempc
'ttempc='tempc
'val=2502.2-2.43089*ttempc'

return('val*1000')

function virtual2(temp,dewp,pres)

*------------------------------------------------- -----------
* Function to return virtual temperature in kelvin given temperature in
* kelvin and dewpoint in kelvin and pressure in mb
*------------------------------------------------- ------------

if (temp > 0 & dewp > 0)
'ppres='pres
'ttemp='temp
'ddewp='dewp
*'vap='satvap2(dewp'-'273.15)
*'mix='mixratio('vap',pres)
*'tempv='virtual(temp,'mix')
'vap='satvap2('ddewp-273.15')
'mix='mixratio('vap','ppres')
'tempv='virtual('ttemp','mix')
else
'tempv=-9999'
endif

return ('tempv')

function virtual(temp,mix)

*------------------------------------------------- -----------
* Function to return virtual temperature given temperature in
* kelvin and mixing ratio in g/g.
*------------------------------------------------- ------------

*'tempv='temp'*(1.0+0.6*'mix')'
'ttemp='temp
'mmix='mix
'tempv=ttemp*(1.0+0.6*mmix)'

return ('tempv')

function satvap2(temp)

*------------------------------------------------- --------------
* Given temp in Celsius, returns saturation vapor pressure in mb
*------------------------------------------------- --------------

*'es=6.112*exp(17.67*'temp'/('temp'+243.5))'
'ttemp='temp
'es=6.112*exp(17.67*ttemp/(ttemp+243.5))'

return('es')

function mixratio(e,p)

*------------------------------------------------- -----
* Given vapor pressure and pressure, return mixing ratio
*-------------------------------------------------------

*'mix=0.622*'e'/('p'-'e')'
'ee='e
'pp='p
'mix=0.622*ee/(pp-ee)'

return('mix')



bullterrier
Участник
Письмо
Пермь
# Дата: 22 Ноя 2012 20:36




Cumulonimbus incus
Участник
Письмо
Кишинёв, Молдова
# Дата: 2 Янв 2013 11:47 - Поправил: Cumulonimbus incus


Совет по отображению новых государственных границ в ГрАДСе.

1. Скачиваем файл http://www.diva-gis.org/data/misc/countries_shp.zi p и распаковываем его в рабочую папку ГрАДСа. По умолчанию путь такой: c:\OpenGrADS\Contents\Resources\SampleDatasets\.

2. Открываем любой файл с данными. После выставления широты и долготы будущей карты вносим в скрипт команду, отключающую вывод основной ГрАДСовой карты:
'set mpdraw off'

3. После отображения переменной даём команду вывода на экран новых границ:
'draw shp countries'

4. Результат:


Пример скрипта:

'reinit'
prompt 'Enter GFS model run date (yyyymmdd) --> '
pull date
prompt 'Enter GFS model run hour (hh) --> '
pull hour
'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_h d'date'/gfs_hd_'hour'z'

'set lat 30 70'
'set lon -20 70'
'set mpdraw off'

'set gxout shade2b'
'd tmp2m-273.15'
'cbarn'
'draw shp countries'


Mesocyclon
Участник
Письмо
Троицк, Челябинск, Челябинская область
# Дата: 21 Янв 2013 12:52


Cumulonimbus incus
Спасибо! Всё работает :)


Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 6 Фев 2013 06:25


Пособие по работе с ГРАДС (на английском) на примере использования японской глобальной модели JMA GSM

http://www.slideshare.net/JMA_447/jma-hr-gsmdatagr ads

Cumulonimbus incus
Участник
Письмо
Кишинёв, Молдова
# Дата: 21 Фев 2013 17:38


Corvus

Хотелось бы у Вас спросить: какие единицы измерения получаются при расчёте лапласиана в ГрАДСе по функции http://www.met.wau.nl/education/atd/practical/func tions/laplace.gs? Мне нужен их перевод в гПа/(300 км)^2 для расчётов параметров конвекции.

Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 21 Фев 2013 17:58


Cumulonimbus incus

Предполагаю, что размерность зависит от шага сетки модели. То есть, например, для ГФС (шаг 0.5 градуса = примерно 50 км) размерность рассчитанного лапласиана будет гПа/(50 км)^2.

Cumulonimbus incus
Участник
Письмо
Кишинёв, Молдова
# Дата: 21 Фев 2013 18:05


Corvus

Маловатые значения получаются - в районе 1,5 10^-8. Буду дальше разбираться.

Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 22 Фев 2013 05:54


Маловатые значения получаются - в районе 1,5 10^-8
Cumulonimbus incus

Возможно, что размерность вообще гПа/км^2.

Судя по этому:

* Horizontal increments
*
"define dy = cdiff(lat,y)*degtorad"
"define dx = cdiff(lon,x)*degtorad"
*
* Define first derivatives
*
"define dhy = cdiff("field",y)"
"define dhx = cdiff("field",x)"
"define dhdy = tcos*dhy/dy"
"define dhdx = dhx/dx"
*
* Define second derivatives
*
"define c2hy = cdiff(dhdy,y)"
"define c2hx = cdiff(dhdx,x)"
"define d2hy = tcos*c2hy/dy"
"define d2hx = c2hx/dx"
*
* Laplacian in spherical coordinates
*
"define laplac = (1/(earthradius*earthradius*tcos*tcos))*(d2hx+d2hy )"

Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 22 Фев 2013 11:12


ГРАДС 2 позволяет делать карты для Гугл.Земля. Надо разбираться. Вот кусочки примеров из разных источников.

To prepare the KML output from GRADS-2 I started by creating a dummy KML file using:

* Zoom in on Nova Scotia
'SET LAT 42.6 47.7'
'SET LON -68.0 -60.9'
* Create Google Earth Dummy Kml file
if (count = 0)
'clear'
'set kml base_kml'
'set gxout kml'
'd mag(UGRD10m,VGRD10m) * 3.6'
endif
* Prep the display for Google Earth map Projections
'clear'
'set parea 0 11 0 8.5'
'set grid off'
'set mproj scaled'
'set mpdraw off '
'set grads off'

and then used the printim command to write out a high resolution GIF or PNG file with the same map projection. This allowed accurate weather forecasts a few days in advance to be calculated for the exact launch site during the planned flight window. Temperature, wind speed and wind direction were plotted as a sequence of PNG or GIF image overlays with the Google Earth KML timespan animation attribute.


Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 22 Фев 2013 11:15


'printim wibble.gif white x720 y381'

From the printim documentation:
If the filename ends with ".png" or ".PNG" then GrADS will
automatically create the image in PNG format
If the filename ends with ".gif" or ".GIF" then GrADS will
automatically create the image in GIF format
If the filename ends with ".jpg" or ".JPG" then GrADS will
automatically create the image in JPEG format
Remember that the printim image has to have the same aspect ratio as
the data grid for it to display properly with the KML file.

>>
>
> how would one make the wibble.kml file?
'set kml wibble'
'set gxout kml'
'd variable'

This gives you two files: wibble.tif and wibble.kml.
You can then edit wibble.kml and change the name of the image file in
the <href> tag to whatever you created with printim.


Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 22 Фев 2013 11:15


The image files created with 'gxout geotiff' and 'gxout kml' are
roughly equivalent to the display inside the plot area of 'gxout
grfill' -- a grid of pixels without anything else. So, if you issue
the following commands:
'set parea 0 11 0 8.5'
'set grid off'
'set mproj scaled'
'set x 0.5 720.5'
'set y 0.5 381.5'
'set mpdraw off'
'set grads off'
'd variable'
'printim x720 y381'

Then you would create a PNG that is roughly equivalent to the TIFF and
you could substitute that file name in the KML and it would work just
as well. In some way, this technique is more flexible because you can
use shaded contours, overlay vectors, etc. and still draw them on
Google Earth.

The TIFF file created with 'gxout kml' has the same geolocation
metadata embedded in it the way the 'gxout geotiff' output does. The
difference is that the 'gxout geotiff' files have floating-point data
values for each pixel and the 'gxout kml' files have color-numbered
index values. I assume there are some GIS applications that would only
need the image and not the data, so that's why I set up a way to
create both. The KML provides a handy way to look at the images with
Google Earth.

For testing, I have been using a free program called Quantum GIS
(qgis.org). Transparency is pretty easy to control within this
application -- you can set any number of pixel values to be
transparent (on a scale of 0-100%). I think Google Earth has some
controls for transparency too. So far, I haven't seen a need to set up
transparency in the GrADS geotiff output. And it's definitely not
clear to me that if I put a transparency mask in the TIFF file, the
GIS applications will know what to do with it. But I'm still a GIS
neophyte...

Please note that some problems with the 'gxout geotiff' output have
been reported to me. I am doing some re-coding and testing and may
have some patches and/or an updated version soon.

Jennifer

p.s. Hooray for TIFF for keeping the standard consistent for all these
years. We need more standards like that.


Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 22 Фев 2013 11:17


This is the technique I have worked out for high quality KML ouput. I
export a dummy kml file + small tiff using gxout kml just for the
Google Earth GroundOverlay <LatLonBox> coordinates. Then I use printim
command to export a larger image in a format like TIFF or PNG.

Then I replace the reference in the KML file from the low resolution
tiff to the new image on the line:
<href>base_kml.tif</href>

Also in my experience Google Earth doesn't respect the transparency in
GrADS TIFF or GIF exports but does in PNG exports.

You can change the image format of the high resolution output to
either PNG, TIFF, etc.. by changing the file extension in the command:
printim full_rez.png x2048 y1024 -t 0

Attached are two Grads Scripts. The script kml_plot_with_countries.gs
will plot political boundary outlines in the exported image. The
script kml_plot_no_outlines.gs disables political boundary outlines
with the grads command set mpdraw off.


Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 22 Фев 2013 11:18


ga-> set gxout kml
ga-> set kml test1
KML output file names:
test1.tif (TIFF image)
test1.kml (KML text file)
KML output type: image
ga-> d tmp2m
Created TIFF image file test1.tif
and complementary KML file test1.kml
ga-> set kml -ln test2
KML output file name:
test2.kml (KML text file)
KML output type: contour lines


Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 22 Фев 2013 11:21


Инфа в мануале

http://www.iges.org/grads/gadoc/gradcomdsetkml.htm l

Cumulonimbus incus
Участник
Письмо
Кишинёв, Молдова
# Дата: 23 Фев 2013 16:38


Corvus

Размерность лапласиана, скорее всего, Па/м^2. Ведь



Нам нужно избавиться от расстояния в знаменателе, чтобы получить привычные значения. Для этого и умножаем на 9*10^8. Прилагаю пример карты (значения похожи на правду):



Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 23 Фев 2013 17:01


Размерность лапласиана, скорее всего, Па/м^2
Cumulonimbus incus

Да, похоже на то :)
Неправ я был с гПа и км.

Cumulonimbus incus
Участник
Письмо
Кишинёв, Молдова
# Дата: 28 Фев 2013 19:10


Corvus

Сейчас как раз делаю новые карты для улучшения прогноза гроз и прочих конв. явлений. Вот первые готовые образцы:



Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 1 Мар 2013 05:28


Cumulonimbus incus

Красивые карты!
Единственное, я бы на второй карте высокие значения удельной влажности раскрасил бы жёлтым и красным цветом, а малые синим. Как на Вашей карте лапласиана.

А то сейчас низкие значения зелёным, высокие синим - это не очень логично. Обычно синий цвет применяют для низких значений.

spralex
Участник
Письмо
г. Конотоп, Сумская обл., Украина
# Дата: 1 Мар 2013 10:45


Cumulonimbus incus

Отличные карты получаются!

TornadoF5
Участник
Письмо
Харьков, Украина. (Игорь)
# Дата: 1 Мар 2013 20:25


Cumulonimbus incus
Карты очень красивые! Спасибо. А что значит параметр "Линии тока в слое 0-30 гПа над поверхностью"? Т.е. это направление и скорость ветра на высоте 10 метров?

Corvus
Автор сайта
Письмо
Владимир (г. Байконур)
# Дата: 1 Мар 2013 20:48


А что значит параметр "Линии тока в слое 0-30 гПа над поверхностью"? Т.е. это направление и скорость ветра на высоте 10 метров?
TornadoF5

Нет, это именно средний ветер в слое 0-30 гПа, то есть примерно в слое 0-300 м над земной поверхностью.

ugrd30_0mb ** 30-0 mb above ground u-component of wind [m/s]
vgrd30_0mb ** 30-0 mb above ground v-component of wind [m/s]

http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_h d20130301/gfs_hd_12z.info

TornadoF5
Участник
Письмо
Харьков, Украина. (Игорь)
# Дата: 1 Мар 2013 20:59


Corvus
Ага, теперь понятно, спасибо.

Cumulonimbus incus
Участник
Письмо
Кишинёв, Молдова
# Дата: 2 Мар 2013 00:16


Corvus, spralex, TornadoF5

Спасибо!

Единственное, я бы на второй карте высокие значения удельной влажности раскрасил бы жёлтым и красным цветом, а малые синим.
Corvus

Это уже, как мне кажется, дело вкуса. У меня прочная ассоциация светлых тонов с сухим воздухом, а тёмных с влажным.

Сейчас сделал пробную карту САРЕ, совмещённую с давлением на уровне моря и сдвигом ветра в слое 0-6 км. Жду комментариев.



Cumulonimbus incus
Участник
Письмо
Кишинёв, Молдова
# Дата: 5 Апр 2013 20:18 - Поправил: Cumulonimbus incus


Прошу помощи у участников Метеоклуба со скриптом для SWEAT. Для своих вычислений я взял формулу внизу. Условные обозначения:

Td850 - точка росы на 850 гПа (град. Цельсия);
TT - индекс Total-Totals;
V850 и V500 - скорости ветра на 850 и 500 гПа (узлы);
dd850 и dd500 - направления ветра на тех же уровнях (градусы).

SWEAT = 12xTd850 + 20x(TT-49) +2x1.94xV850 + 1.94xV500 + 125x(S+0.2), где S = sin(dd500-dd850). Но этот синус считается только при следующих условиях:

1) 130°=<dd850=<250°
2) 210°=<dd500=<310°
3) dd500-dd850 > 0°
4) V850 and V500 >= 15 узлов

Исходя их этого, составил такой код:

function sweat(td850,tti,v850,v500,dd850,dd500)
'ptd850 = maskout(12*'td850',12*'td850')'
'ptd850 = const(ptd850,0,-u)'
'ptti = maskout(20*('tti'-49),20*('tti'-49))'
'ptti = const(ptti,0,-u)'
'pv850 = 2*'v850''
'pv500 = 'v500''
'deltadd = 0.0174532925*('dd500'-'dd850')'
'deltadd = maskout(deltadd,'dd850'-130)'
'deltadd = maskout(deltadd,250-'dd850')'
'deltadd = maskout(deltadd,'dd500'-210)'
'deltadd = maskout(deltadd,310-'dd500')'
'deltadd = maskout(deltadd,deltadd)'
'deltadd = maskout(deltadd,'v850'-15)'
'deltadd = maskout(deltadd,'v500'-15)'
'ps = 125*(sin(deltadd)+0.2)'
'ps = const(ps,0,-u)'
'sweat = ptd850+ptti+pv850+pv500+ps'
return


Также воспользовался готовым кодом с библиотеки ГрАДСа:

ftp://cola.gmu.edu/grads/scripts/sweat_index.gs

При расчётах появились существенные расхождения моего кода с библиотечным в пятом слагаемом (с синусом). Исходные данные следующие:

dd850 = 208.559 град.;
dd500 = 252.755 град.;
v850 = 13.7435 узлов;
v500 = 37.8981 узлов.

Так как v850 < 15 узлов, то весь член с синусом обнуляется, что и выдаёт мой скрипт. А чужой код рассчитывает 112.139, будто такого условия и нет. Кто прав?

TornadoF5
Участник
Письмо
Харьков, Украина. (Игорь)
# Дата: 5 Апр 2013 21:28


Cumulonimbus incus
Я в математике и программировании не силён, поэтому к сожалению в этой ситуации помочь ничем не могу... :(

<< . 1 . 2 . 3 . 4 . 5 . 6 . 7 . 8 . 9 . 10 . >>
Ваш ответ

          Отменить *Что это?

 » Логин  » Пароль 
 
 
Полезная информация:



Поддержка: miniBB forum software © 2001-2017