mini_ker.texi 141 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199320032013202320332043205320632073208320932103211321232133214321532163217321832193220322132223223322432253226322732283229323032313232323332343235323632373238323932403241324232433244324532463247324832493250325132523253325432553256325732583259326032613262326332643265326632673268326932703271327232733274327532763277327832793280328132823283328432853286328732883289329032913292329332943295329632973298329933003301330233033304330533063307330833093310331133123313331433153316331733183319332033213322332333243325332633273328332933303331333233333334333533363337333833393340334133423343334433453346334733483349335033513352335333543355335633573358335933603361336233633364336533663367336833693370337133723373337433753376337733783379338033813382338333843385338633873388338933903391339233933394339533963397339833993400340134023403340434053406340734083409341034113412341334143415341634173418341934203421342234233424342534263427342834293430343134323433343434353436343734383439344034413442344334443445344634473448344934503451345234533454345534563457345834593460346134623463346434653466346734683469347034713472347334743475347634773478347934803481348234833484348534863487348834893490349134923493349434953496349734983499350035013502350335043505350635073508350935103511351235133514351535163517351835193520352135223523352435253526352735283529353035313532353335343535353635373538353935403541354235433544354535463547354835493550355135523553355435553556355735583559356035613562356335643565356635673568356935703571357235733574357535763577357835793580358135823583358435853586358735883589359035913592359335943595359635973598359936003601360236033604360536063607360836093610361136123613361436153616361736183619362036213622362336243625362636273628362936303631363236333634363536363637363836393640364136423643364436453646364736483649365036513652365336543655365636573658365936603661366236633664366536663667366836693670367136723673367436753676367736783679368036813682368336843685368636873688368936903691369236933694369536963697369836993700370137023703370437053706370737083709371037113712371337143715371637173718371937203721372237233724372537263727372837293730373137323733373437353736373737383739374037413742374337443745374637473748374937503751375237533754375537563757375837593760376137623763376437653766376737683769377037713772377337743775377637773778377937803781378237833784378537863787378837893790379137923793379437953796379737983799380038013802380338043805380638073808380938103811381238133814381538163817381838193820382138223823382438253826382738283829383038313832383338343835383638373838383938403841384238433844384538463847384838493850385138523853385438553856385738583859386038613862386338643865386638673868386938703871387238733874387538763877387838793880388138823883388438853886388738883889389038913892389338943895389638973898389939003901390239033904390539063907390839093910391139123913391439153916391739183919392039213922392339243925392639273928392939303931393239333934393539363937393839393940394139423943394439453946394739483949395039513952395339543955395639573958395939603961396239633964396539663967396839693970397139723973397439753976397739783979398039813982398339843985398639873988398939903991399239933994399539963997399839994000400140024003400440054006400740084009401040114012401340144015401640174018401940204021402240234024402540264027402840294030403140324033403440354036403740384039404040414042404340444045404640474048404940504051405240534054405540564057405840594060406140624063406440654066406740684069407040714072407340744075407640774078407940804081408240834084408540864087408840894090409140924093409440954096409740984099410041014102410341044105410641074108410941104111411241134114411541164117411841194120412141224123412441254126412741284129413041314132413341344135413641374138413941404141414241434144414541464147414841494150415141524153415441554156415741584159416041614162416341644165416641674168416941704171417241734174417541764177417841794180418141824183418441854186418741884189419041914192419341944195419641974198419942004201420242034204420542064207420842094210421142124213421442154216421742184219
  1. \input texinfo @c -*-texinfo-*-
  2. @ignore
  3. verifier si on a bien les nouveautes de la 1.02
  4. 1/ sensy_to_var
  5. 2/ Borel et Boreleig OK
  6. 3/ v.p. param OK
  7. 4/ Helps et listing divers OK
  8. 5/ Revoir Probes avec Pat => et code !
  9. 6/ data doivent-elles etre specifique de kalman?
  10. 7/ symbolic en premier dans Borel a faire
  11. @end ignore
  12. @setfilename mini_ker.info
  13. @include version.texi
  14. @c @set myversion @value{VERSION}
  15. @set myversion 102
  16. @set myurl @url{http://www.environnement.ens.fr/@/perso/@/dumas/@/mini_ker/@/software.html}
  17. @macro Minik{}
  18. Miniker
  19. @end macro
  20. @settitle @Minik{} @value{myversion} manual
  21. @syncodeindex fn vr
  22. @dircategory Miscellaneous
  23. @direntry
  24. * Miniker: (mini_ker). The mini_ker modeling tool.
  25. @end direntry
  26. @copying
  27. @noindent Copyright (C) 2004, 2005, 2006, 2007 Alain Lahellec@*
  28. Copyright (C) 2004, 2005, 2006, 2007 Patrice Dumas@*
  29. Copyright (C) 2004, St@'ephane Hallegatte
  30. @quotation
  31. Permission is granted to copy, distribute and/or modify this document
  32. under the terms of the GNU Free Documentation License, Version 1.1 or
  33. any later version published by the Free Software Foundation; with no
  34. Invariant Sections, with no Front-Cover text and with no Back-Cover Text.
  35. A copy of the license is included in the section entitled ``GNU Free
  36. Documentation License.''
  37. @end quotation
  38. @end copying
  39. @titlepage
  40. @title @Minik{} manual
  41. @subtitle for @Minik{} version @value{myversion}, @value{UPDATED}
  42. @author The TEF Collaboration
  43. @c @author Alain Lahellec
  44. @c @author Patrice Dumas
  45. @c @author St@'ephane Hallegatte
  46. @page
  47. @vskip 0pt plus 1filll
  48. @insertcopying
  49. @end titlepage
  50. @ifset texi2html
  51. @node Top
  52. @top @Minik{} @value{myversion} manual
  53. @c @insertcopying
  54. @end ifset
  55. @ifnottex
  56. @node Top
  57. @top @Minik{} @value{myversion} manual
  58. @strong{By: The TEF Collaboration}
  59. @insertcopying
  60. @end ifnottex
  61. @menu
  62. * Introduction::
  63. * TEF overview::
  64. * A model with @Minik{}::
  65. * Advanced programming::
  66. * Dynamic system analysis::
  67. * Advanced use of @Minik{} with make::
  68. Indices
  69. * Concepts index::
  70. * Variables macros and functions index::
  71. Appendices
  72. * Installation::
  73. * Cmz directives reference::
  74. @ignore
  75. * Resolution method::
  76. * @Minik{} macros::
  77. @end ignore
  78. * Copying This Manual:: The GNU Free Documentation License.
  79. @end menu
  80. @contents
  81. @node Introduction
  82. @unnumbered Introduction
  83. @cindex TEF
  84. @cindex cells
  85. @cindex transfers
  86. @cindex ZOOM
  87. @cindex mortran
  88. @Minik{} is a modeling tool, built especially in order to implement
  89. models written following the @acronym{TEF,Transfer Evolution Formalism}
  90. formalism, a mathematical framework for system analysis and
  91. simulation. @Minik{} allows for timewise simulation, system analysis,
  92. adjoint computation, Kalman filtering and more.
  93. @Minik{} uses a fortran preprocessor, @command{mortran}, designed in the
  94. 1970's, to ease model writing using dedicated specific languages.
  95. For example partial derivatives are
  96. symbolicaly determined by @command{mortran} macros in @Minik{}.
  97. For the selection of
  98. another compile-time features, another set of preprocessor directives,
  99. the @dfn{cmz directives}, are used. In most cases the user does not need to
  100. know anything about that preprocessing that occurs behind the scene,
  101. he simply writes down the equations of his model and he is done.
  102. @c All partial derivatives needed to solve the TEF system are automatically
  103. @c determined during the pre-compilation stage.
  104. @c Once all models written down and initial conditions
  105. @c given using a pseudo-Fortran type of language, the model is ready to run.
  106. @c The language developed to get automatic symbolic partial derivatives
  107. @c uses the Fortran pre-compiler @command{mortran}, designed in the 1970's.
  108. A comprehensive description
  109. of the @acronym{TEF} formalism in available on
  110. @url{http://www.lmd.jussieu.fr/ZOOM/doc/tef-GB-partA5.pdf}).
  111. The @Minik{} software is a reduced version of
  112. @uref{http://www.lmd.jussieu.fr@//zoom,@strong{ZOOM}},
  113. that can only handle a hundreds of variables, but is much easier to use.
  114. @menu
  115. * Intended audience::
  116. * Reading guide::
  117. * Other Manuals and documentation::
  118. @end menu
  119. @node Intended audience
  120. @unnumberedsec Intended audience
  121. The reader should have notions in system dynamics.
  122. @c and understand the basis of the TEF.
  123. Moreover a minimal knowledge of programmation and fortran is
  124. required. What is required is a basic understanding of variable types,
  125. affectation and fortran expressions.
  126. @node Reading guide
  127. @unnumberedsec Reading guide
  128. The first chapter is a brief overview of the @acronym{TEF}.
  129. The following describes how to write, compile and run a model in @Minik{}
  130. in its basic and comprehensive syntax.
  131. @c Reading the sections of this chapter up to the section
  132. @c @emph{Symbolic model description} is required to know the
  133. @c syntax of model description in @Minik{}.
  134. Reading up to the section
  135. @emph{Controlling the run} is required to be able to use @Minik{}.
  136. In this section it is assumed that @Minik{} is properly setup. The
  137. installation instructions are in the appendix at
  138. @ref{Installation}.
  139. @c 2 programming environment to compile the model are available, with cmz
  140. @c and make, you can skip the method description you are not interested in.
  141. @c A reference for the usefull cmz directives is also in the appendix
  142. @c (@pxref{Cmz directives reference}).
  143. @c You should also
  144. @c read the following section, @emph{Symbolic model description} which presents an
  145. @c alternate syntax for model description, such that you can choose what you
  146. @c prefer.
  147. The next chapter describes advanced features, first a general introduction to
  148. features settings and then a description of other model description related
  149. features.
  150. The next chapter describes system analysis tools available with @Minik{}. The
  151. sections are independant and each describes how to use a specific feature. If
  152. you plan on using these features, you should also read
  153. @ref{Selecting features, , Overview of feature setting}.
  154. A final chapter describes advanced features in a development environment
  155. using make,
  156. In the appendix the instructions for the installation are described
  157. (@pxref{Installation}).
  158. @node Other Manuals and documentation
  159. @unnumberedsec Other Manuals and documentation
  160. A programmers'Manual is available (in French), and can be asked for to
  161. any member of the collabration. See additional documents in
  162. @url{http://www.lmd.jussieu.fr/Zoom/doc} or ask for Research
  163. texts and articles to members.
  164. @node TEF overview
  165. @chapter An overview of the @acronym{TEF} formalism
  166. The @acronym{TEF, Transfer Evolution Formalism} is based on partitionning
  167. and recoupling of model subsystems. It allows the study of the coupling
  168. between subsystems by the means of linearization and time discretization.
  169. @menu
  170. * Cell and Transfer::
  171. * Linearization and discretization::
  172. @end menu
  173. @node Cell and Transfer
  174. @section Cell and Transfer equations
  175. In the @acronym{TEF}, a model is
  176. mathematically represented by a set of equations corresponding
  177. to two kinds objects:
  178. @enumerate
  179. @item Cells which are elementary models and correspond to evolution equations
  180. such as:
  181. @ifset latex
  182. @tex
  183. \begin{eqnarray}
  184. \partial_t \eta (t) &=& g(\eta(t),\varphi(t))
  185. \label{cells}
  186. \end{eqnarray}
  187. @end tex
  188. @end ifset
  189. @ifclear latex
  190. @tex
  191. $$\partial_t \eta (t) = g(\eta(t),\varphi(t))$$
  192. @end tex
  193. @ifnottex
  194. @noindent @math{d eta(t)/d t = g(eta(t),phi(t))}
  195. @end ifnottex
  196. @end ifclear
  197. Vector @math{\eta} represent the state variables of cells and
  198. the vector @math{\varphi} represent the dependent
  199. boundary conditions, @i{i.e.} the
  200. variables considered as boundary conditions by a cell, but depending upon
  201. the complete model state. This dependent boundary conditions are
  202. required to make the cells correspond to well-posed problems.
  203. @c FIXME acceptable?
  204. These variables are often called state variables, and prognostic
  205. variables in meteorology.
  206. @item Transfers which are determined by constraint equations such as:
  207. @ifset latex
  208. @tex
  209. \begin{eqnarray}
  210. \varphi(t) &=& f(\eta(t),\varphi(t))
  211. \label{transfer}
  212. \end{eqnarray}
  213. @end tex
  214. @end ifset
  215. @ifclear latex
  216. @tex
  217. $$
  218. \varphi(t) = f(\eta(t),\varphi(t))
  219. $$
  220. @end tex
  221. @ifnottex
  222. @noindent @math{phi(t) = f(eta(t),phi(t))}
  223. @end ifnottex
  224. @end ifclear
  225. These equations are often called algebraic equations, and in meteorology
  226. diagnostic equations.
  227. @end enumerate
  228. @node Linearization and discretization
  229. @section Linearization and discretization in the @acronym{TEF}
  230. The relations between sub-systems is excessively difficult to exhibit when
  231. having to cope with non-linear system. In the @acronym{TEF}, the
  232. @acronym{TLS, Tangent Linear System} is constructed along the trajectory.
  233. One considers the system over a small portion along the trajectory, say
  234. between @math{t} and @math{t + \delta t}. The variation @math{\delta \eta}
  235. of @math{\eta} and @math{\delta \varphi} of @math{\varphi} is obtained
  236. through a Pad@'e approximation of the state-transition matrix. The final
  237. form of the algebraic system is closed to the classical Crank-Nicolson scheme:
  238. @c FIXME PAd'e? od Taylor?
  239. @c through a Taylor expansion followed by time integration.
  240. @c A time scheme is then applied to the @acronym{TLS} (a Crank-Nicholson scheme),
  241. @c to obtain an algebraic system describing the relationships between
  242. @c variations of transfers and cells variables:
  243. @ifset latex
  244. @tex
  245. $$\left(\begin{array}{cc}
  246. A & B\\
  247. -C^+ & I-D
  248. \end{array}\right) \left(\begin{array}{c}
  249. \delta \eta\\
  250. \delta \varphi
  251. \end{array}\right) = \left(\begin{array}{c}
  252. \Gamma\\
  253. \Omega \end{array}\right)$$
  254. @end tex
  255. @end ifset
  256. @ifclear latex
  257. @tex
  258. $$\pmatrix{A & B\cr
  259. -C^+ & I-D\cr} \pmatrix{\delta \eta\cr
  260. \delta \varphi\cr} = \pmatrix{\Gamma\cr
  261. \Omega\cr}$$
  262. @end tex
  263. @end ifclear
  264. The blocks appearing in the Jacobian matrix are constructed with partial derivative
  265. of @math{f} and @math{g}, and with @math{\delta t}. From this system the
  266. elimination of @math{\delta \eta} leads to another formulation giving
  267. the coupling between transfers, and allows for the @math{\delta \varphi}
  268. computation. The @math{\delta \varphi} value is then substitued in
  269. @math{\delta \eta} to complete the time-step solving process.
  270. @node A model with @Minik{}
  271. @chapter @Minik{} model programming
  272. @cindex sequences
  273. @Minik{} works by combining the model specification code given by
  274. the user and other source files provided in the package. The code is
  275. assembled, preprocessed, compiled, linked and the resulting program
  276. can be run to produce the model trajectory and dynamic analysis.
  277. The code provided in the package contains a principal program, some usefull
  278. subroutines and pieces of code called @dfn{sequences} combined with the
  279. different codes. Among these sequences some hold the code describing the model
  280. and are to be written by the user (sequences are similar to
  281. Fortran include files).
  282. @menu
  283. * Structure of the code::
  284. * A model description::
  285. * Setting and running a model::
  286. * Controlling the run::
  287. @end menu
  288. @node Structure of the code
  289. @section General structure of the code
  290. @cindex sequence
  291. @cindex zinit, general
  292. The sequences used to enter model description hold the @c vector dimensions,
  293. mathematical formulae for each cell and transfer component, dedicated
  294. derived computations, and time-step
  295. steering. During the code generation stage,
  296. cmz directives are preprocessed, then the user pseudo-Fortran
  297. instructions are translated by @command{mortran} using macros designed to
  298. generate in particular all Fortran instructions that compute the Jacobian
  299. matrices used in @acronym{TEF} modelling.
  300. @c A first users' sequence to program is: @file{dimetaphi} where the model
  301. @c dimensions are given, for the two vector-array @code{eta(.)} for cells
  302. @c and @code{ff(.)} for transfers (@pxref{Model size,,Entering model size}).
  303. The sequence @file{zinit} contains the mathematical formulation of the model
  304. (@pxref{Model equation and parameters, Entering model equation and parameters}).
  305. Another sequence, @file{zsteer}, is merged at
  306. the end of the time step advance of the simulation, where the user can
  307. monitor the time step values and printing levels, and perform particular
  308. computations etc.
  309. (@pxref{End of time step, ,Executing code at the end of each time step}).
  310. @node A model description
  311. @section @Minik{} programming illustrated
  312. @cindex TEF
  313. The general @acronym{TEF} system writes:
  314. @ifset latex
  315. @tex
  316. \begin{eqnarray}
  317. \partial_t \eta (t) &=& g(\eta(t),\varphi(t)) \nonumber\\
  318. \varphi(t) &=& f(\eta(t),\varphi(t))
  319. \label{tef}
  320. \end{eqnarray}
  321. @end tex
  322. @end ifset
  323. @ifclear latex
  324. @tex
  325. $$\eqalign{\partial_t \eta (t) &= g(\eta(t),\varphi(t))\cr
  326. \varphi(t) &= f(\eta(t),\varphi(t))\cr
  327. }$$
  328. @end tex
  329. @ifnottex
  330. @noindent @math{d eta(t)/d t = g(eta(t),phi(t))@*
  331. phi(t) = f(eta(t),phi(t))}
  332. @end ifnottex
  333. @end ifclear
  334. To illustrate the model description in @Minik{} a simple predator-prey
  335. model of Lotka-Volterra is used.
  336. This model can be written in the following @acronym{TEF} form:
  337. @ifset latex
  338. @tex
  339. \begin{equation}
  340. \left\{
  341. \begin{array}{cc}
  342. \partial_t \eta _{prey} =& a \eta _{prey} - a \varphi _{meet} \\
  343. \partial_t \eta _{pred} =& -c \eta _{pred} + c \varphi _{meet} \nonumber\\
  344. \end{array}
  345. \right.
  346. \end{equation}
  347. \begin{equation}
  348. \varphi _{meet} = \eta _{prey}\eta _{pred}
  349. \label{pred}
  350. \end{equation}
  351. @end tex
  352. @end ifset
  353. @ifclear latex
  354. @tex
  355. $$\left\{\eqalign{\partial_t \eta _{prey} &= a \eta _{prey} - a \varphi _{meet} \cr
  356. \partial_t \eta _{pred} &= -c \eta _{pred} + c \varphi _{meet}\cr}\right.$$
  357. @end tex
  358. @tex
  359. $$\varphi _{meet} = \eta _{prey}\eta _{pred}$$
  360. @end tex
  361. @ifnottex
  362. @noindent @math{d eta_prey(t)/d t = a * eta_prey - a * phi_meet@*
  363. d eta_pred(t)/d t = -c * eta_pred +c * phi_meet}
  364. @noindent @math{phi_meet = eta_prey * eta_pred}
  365. @end ifnottex
  366. @end ifclear
  367. with two cell equations, @i{i.e}. state evolution of the prey and predator
  368. groups, and one transfer accounting for the meeting of individuals of
  369. different group.
  370. @menu
  371. * A note about mortran and cmz directives::
  372. * Model equation and parameters::
  373. @end menu
  374. @node A note about mortran and cmz directives
  375. @subsection All you need to know about mortran and cmz directives
  376. @cindex mortran
  377. The first stage of code generation consists in cmz directives preprocessing.
  378. Cmz directives are used for conditional selection of features, and sequence
  379. inclusion. At that point you don't need to know anything about these
  380. directives. They are only usefull if you want to take advantage of advanced
  381. features
  382. (@pxref{Programming with cmz directives}).
  383. The code in sequences is written in Mortran and the second stage of code
  384. generation consists in mortran macro expansion. The mortran language is
  385. described
  386. in its own manual, here we only explain the very basics which is all you need
  387. to know to use @Minik{}. Mortran basic instructions are almost Fortran,
  388. the differences are the following:
  389. @itemize @bullet
  390. @item The code is free-form, and each statement should end with a semi-colon
  391. @code{;}.
  392. @item Comments may be introduced by an exclamation mark @code{!} at the
  393. beginning of a line, or appear within double quotes @code{"} in a single line.
  394. @item It is possible to use blocs, for @code{do} or @code{if} statement
  395. for example, and they are enclosed within brackets @samp{<} and @samp{>}.
  396. To be in the safe side, a semi-colon @code{;} should be added after a
  397. closng bracket @code{>}.
  398. @end itemize
  399. The following fictious code is legal mortran:
  400. @example
  401. real
  402. param;
  403. param = 3.; ff(1) = ff(3)**eta(1); "a comment"
  404. ! a line comment
  405. do inode=1,n_node <eta_move(inode)=0.01; eta_speed(inode)=0.0;>;
  406. @end example
  407. Thanks to mortran the model code is very simply specified, as you'll
  408. see next.
  409. @node Model equation and parameters
  410. @subsection Entering model equation and parameters
  411. @cindex @file{zinit}
  412. @vindex dt
  413. @vindex time
  414. @vindex nstep
  415. @vindex modzprint
  416. The model equation and parameters and some @Minik{} parameters are entered in
  417. the @file{zinit} sequence. The whole layout of the model is given
  418. before detailing the keywords.
  419. @example
  420. !%%%%%%%%%%%%%%%%%%%%%%
  421. ! Parameters
  422. !%%%%%%%%%%%%%%%%%%%%%%
  423. real apar,bpar; "optional Fortran type declaration"
  424. ! required parameters
  425. dt=.01; "initial time-step"
  426. nstep=10 000; "number of iterations along the trajectory"
  427. time=0.; "time initialisation "
  428. ! model parameters
  429. apar = 1.5;
  430. cpar = 0.7;
  431. ! misceallaneous parameters
  432. modzprint = 1000; "printouts frequency"
  433. print*,'***************************************';
  434. print*,'Lotka-Volterra model with parameters as:';
  435. z_pr: apar,bpar;
  436. print*,'***************************************';
  437. !%%%%%%%%%%%%%%%%%%%%%%
  438. ! Transfer definition
  439. !%%%%%%%%%%%%%%%%%%%%%%
  440. ! rencontre (meeting)
  441. set_Phi
  442. < var: ff_interact, fun: f_interact = eta_prey*eta_pred;
  443. >;
  444. !%%%%%%%%%%%%%%%%%%%%%%
  445. ! Cell definition
  446. !%%%%%%%%%%%%%%%%%%%%%%
  447. set_eta
  448. < var: eta_prey, fun: deta_prey = apar*eta_prey - apar*ff_interact;
  449. var: eta_pred, fun: deta_pred = - cpar*eta_pred + cpar*ff_interact;
  450. >;
  451. !%%%%%%%%%%%%%%%%%%%%%%
  452. ! Initial states
  453. !%%%%%%%%%%%%%%%%%%%%%%
  454. eta_prey = 1.;
  455. eta_pred = 1.;
  456. ;
  457. OPEN(50,FILE='title.tex',STATUS='UNKNOWN'); "title file"
  458. write(50,5000) apar,cpar;
  459. 5000;format('Lotka-Volterra par:',2F4.1);
  460. @end example
  461. @subsubheading Variables and model parameters
  462. The following variables are mandatory:
  463. @table @code
  464. @item dt
  465. The time step.
  466. @item time
  467. Model time initialisation.
  468. @item nstep
  469. Number of iterations along the trajectory.
  470. @end table
  471. There are no other mandatory variables. Some optional variables are used
  472. to monitor the printout and ouput of results of the code.
  473. As an example, the variable @code{modzprint} is used to set
  474. the frequency of the printout of the model matrix and vectors during the
  475. run (@pxref{Controlling the printout and data output}).
  476. User's defined variable and Fortran or Mortran instructions can always be
  477. added for intermediate calculus. To avoid conflict with the variables of the
  478. @Minik{} code, the rule is that a users symbol must not have characters
  479. @samp{o}
  480. in the first two symbol characters.
  481. In the predator-prey example there are two model parameters. The fortran
  482. variables are called here @code{apar} for @math{a} and @code{cpar} for @math{c}.
  483. If a Fortan type definition is needed, it should be set at the very beginning
  484. of @file{zinit}. The predator-prey code variable initializations finally reads
  485. @example
  486. @group
  487. !%%%%%%%%%%%%%%%%%%%%%%
  488. ! Parameters
  489. !%%%%%%%%%%%%%%%%%%%%%%
  490. real apar,bpar; "optional Fortran type declaration"
  491. dt=.01;
  492. nstep=10 000;
  493. time=0.;
  494. ! model parameters
  495. apar = 1.5;
  496. cpar = 0.7;
  497. modzprint = 1000;
  498. @end group
  499. @end example
  500. @subsubheading Model equations
  501. @anchor{Model equations}
  502. @findex set_Phi
  503. @findex set_eta
  504. @vindex var:
  505. @vindex fun:
  506. @vindex eqn:
  507. The model equations for cells and model equations for transferts are
  508. entered in two mortran blocks, one for the transferts, the other for the
  509. cell components. The model equations for cells are entered into a
  510. @code{set_eta} block, and the transfer equations are entered into a
  511. @code{set_phi} block.
  512. In each block the couples variable-function are specified. For
  513. transfers the function defines the transfer itself while for cells
  514. the function describes the cell evolution. The variable is specified
  515. with @code{var:}, the function is defined with @code{fun:}.
  516. In the case of the predator-prey model, the transfer variable
  517. associated with @math{\varphi_{meet}} could be called @code{ff_interact}
  518. and the transfer definition would be given by:
  519. @example
  520. set_Phi
  521. < var: ff_interact, fun: f_interact = eta_prey*eta_pred;
  522. >;
  523. @end example
  524. The two cell equations of the predator-prey model, with name
  525. @code{eta_prey} for the prey (@math{\eta_{prey}}) and @code{eta_pred}
  526. for the predator (@math{\eta_{pred}}) are:
  527. @example
  528. set_eta
  529. < var: eta_prey, fun: deta_prey = apar*eta_prey - apar*ff_interact;
  530. var: eta_pred, fun: deta_pred = - cpar*eta_pred + cpar*ff_interact;
  531. >;
  532. @end example
  533. The @samp{;} at the end of the mortran block is important.
  534. @page
  535. The whole model equations are setup with:
  536. @example
  537. @group
  538. !%%%%%%%%%%%%%%%%%%%%%%
  539. ! Transfer definition
  540. !%%%%%%%%%%%%%%%%%%%%%%
  541. ! rencontre (meeting)
  542. set_Phi
  543. < var: ff_interact, fun: f_interact = eta_prey*eta_pred;
  544. >;
  545. !%%%%%%%%%%%%%%%%%%%%%%
  546. ! Cell definition
  547. !%%%%%%%%%%%%%%%%%%%%%%
  548. set_eta
  549. < var: eta_prey, fun: deta_prey = apar*eta_prey - apar*ff_interact;
  550. var: eta_pred, fun: deta_pred = - cpar*eta_pred + cpar*ff_interact;
  551. >;
  552. @end group
  553. @end example
  554. Whenever the user is not concerned with giving a specific name to a
  555. function, it is possible to specify the equation only with
  556. @code{eqn:}. Therefore the user may replace an instruction as:
  557. @example
  558. var: ff_dump,
  559. fun: f_dump = - rd*(eta_speed - eta_speed_limiting);
  560. @end example
  561. with:
  562. @example
  563. eqn: ff_dump = - rd*(eta_speed - eta_speed_limiting);
  564. @end example
  565. In that case, the unnamed function will take the name of the defined
  566. variable preceded by the @samp{$} sign: @code{$ff_dump}.
  567. @subsubheading Starting points
  568. @cindex starting point
  569. The cells equations require state initial conditions. In some case, the
  570. transfers may also need starting points although they are determined from
  571. the cell values.
  572. In the predator-prey model the starting points for cells are:
  573. @example
  574. ! initial state
  575. ! -------------
  576. eta_prey = 1.;
  577. eta_pred = 1.;
  578. @end example
  579. When there is a non trivial implicit relationship between the transfers
  580. in the model, it may be usefull or even necessary to set some
  581. transfers to non-zero values. This difficulty is only relevant for the very
  582. first step of the simulation and will be used as a
  583. first guess of @math{\varphi}. The uninitialized transfers having
  584. a default compiler-dependant (often zero) value, an initialization
  585. to another value may help avoiding singular functions or matrix and
  586. ensure convergence in the Newton algorithm used to solve the transfer implicit
  587. equation.
  588. @ignore
  589. Indeed a good starting
  590. point for the transfers may help finding their value at the first time step
  591. (to help
  592. avoiding a singular matrix during the research of the first transfers and
  593. ensure convergence), when the implicit equation defining transfers is solved:
  594. @ifset latex
  595. @tex
  596. \begin{eqnarray}
  597. \varphi(t) &=& f(\eta(t),\varphi(t))
  598. \label{transfer}
  599. \end{eqnarray}
  600. @end tex
  601. @end ifset
  602. @ifclear latex
  603. @tex
  604. $$
  605. \varphi(t) = f(\eta(t),\varphi(t))
  606. $$
  607. @end tex
  608. @ifnottex
  609. @noindent @math{phi(t) = f(eta(t),phi(t))}
  610. @end ifnottex
  611. @end ifclear
  612. @end ignore
  613. @subsubheading The cell and transfer arrays
  614. @vindex eta(.)
  615. @vindex ff(.)
  616. @vindex mp
  617. @vindex np
  618. Sometime it is easier to iterate over an array than to use the
  619. cell or transfer variable name. This is possible because there is a
  620. correspondence between the variable names
  621. and the fortran array @code{eta(.)} for the cell variables and
  622. the fortran array @code{ff(.)} for the transfer variables@footnote{In fact
  623. the variables names are transformed into fortran array elements
  624. by mortran generated macros, so the symbolic names defined in the
  625. mortran blocks never appears in the generated fortran code, they are
  626. replaced by the fortran arrays.}.
  627. The index of the variable is determined by the order of appearance in
  628. the variable definition blocks. It is reminded in the output, as
  629. explained later (@pxref{Simulation and output}).
  630. The number of cells is in the integer @code{np} variable, and the
  631. number of transfer is in the integer @code{mp} variable.
  632. @subsubheading title file
  633. @anchor{Title file}
  634. @cindex title file
  635. @cindex @file{title.tex}
  636. For some graphics generation, a file with name @file{title.tex} is required
  637. which sets the title. The following instructions take care of that:
  638. @example
  639. OPEN(50,FILE='title.tex',STATUS='UNKNOWN');
  640. write(50,5000) apar,cpar;
  641. 5000;format('Lotka-Volterra par:',2F4.1);
  642. close(50);
  643. @end example
  644. In that case the parameter values are written down, to differenciate between
  645. different runs. This step is in general not needed.
  646. @c The correspondence with basic components are printed out at execution
  647. @c time as explained in @ref{Simulation and output,,
  648. @c Running a simulation and using the output}. Also, a @file{Model.hlp} is
  649. @c generated that recalls the basic names and equations of the model.
  650. @c It may be noted that whenever
  651. @c the order of variable-functions is the same between indexed declaration and
  652. @c symbolic, the two generated Fortran code are almost identical.
  653. @node Setting and running a model
  654. @section Setting and running a model
  655. In this section it is assumed that a programming environment has been
  656. properly setup. This environment may use either cmz or make to drive
  657. the preprocessing and compilation.
  658. You can skip the part related with the environment you don't intend to use.
  659. For instructions regarding the
  660. installation, see @ref{Installation}.
  661. @menu
  662. * Setting up a model with cmz::
  663. * Setting up a model with make::
  664. * Simulation and output::
  665. * Graphics::
  666. @end menu
  667. @node Setting up a model with cmz
  668. @subsection Setup a model and compile with cmz
  669. @cindex @command{mod}
  670. @cindex @file{$zinit}
  671. @cindex @file{$dimetaphi}
  672. The user defined sequences are @samp{KEEP} in the
  673. cmz world. The most common organization is to have a cmz file in a
  674. subdirectory of the directory containing the @file{mini_ker.cmz}
  675. cmz file. In this
  676. cmz file there should be a @samp{PATCH} called @samp{zinproc}
  677. with the KEEPs within the patch. The KEEP must be called @file{$zinit}.
  678. @c and @file{$dimetaphi}.
  679. From within cmz in the directory of your model the source extraction,
  680. compilation and linking will be triggered by a @command{mod} command. This macro
  681. uses the @file{selseq.kumac} information to find the @file{mini_ker.cmz}
  682. cmz file.
  683. @command{mod}
  684. shall create a directory with the same name than the cmz file,
  685. @file{mymodel/} in our example. In this directory there is another
  686. directory @file{cfs/} containing the sources extracted from the cmz file.
  687. The file @file{mymodel_o.tmp} contains all the mortran code generated
  688. by cmz with the sequences substituted, including the @file{$zinit}. @c and
  689. @c @file{$dimetaphi} sequences (assembled code).
  690. The fortran produced by the preprocessing and
  691. splitting of this file is in files with the traditional @samp{.f} suffix.
  692. The principal program is in @file{principal.f}. An efficient way of getting
  693. familiar with mini_ker methods is looking at the @file{mymodel_o.tmp} where
  694. all sequences and main Mortran instructions are gathered. Symbolic derivation
  695. @c FIXME pas ici la symbolic derivation
  696. is noted as @code{F_D(expression)(/variable)}, and the resulting Fortran code
  697. is in @file{principal.f}.
  698. @command{mod} also triggers compilation and linking. The object files are in
  699. the same @file{cfs/} directory and the executable is in the @file{mymodel/}
  700. directory, with name @file{mymodel.exe}.
  701. @node Setting up a model with make
  702. @subsection Setup a model and compile with make
  703. @cindex compilation
  704. @c @cindex @file{dimetaphi.mti}
  705. @cindex @file{zinit.mti}
  706. @vindex model_file_name
  707. With make, the sequences are files ending with @samp{.mti} (for
  708. mortran include files),
  709. called, for example, @file{zinit.mti}.
  710. @c and @file{dimetaphi.mti}.
  711. They are included by
  712. @command{mortran} in other source files. You also need a @file{Makefile}
  713. to drive the compilation of the model.
  714. If you don't need additional code or libraries to be linked with your
  715. model you have two alternatives.
  716. @enumerate
  717. @item
  718. The simplest alternative is to run
  719. the @command{start_miniker} script with the model file name as argument.
  720. It should copy a @file{zinit.mti} file
  721. ready to be edited and a Makefile ready to compile the model. For
  722. the predator prey model, for example, you could run
  723. @example
  724. $ start_miniker predator
  725. @end example
  726. @item
  727. Otherwise you can copy the Makefile from @file{template/Makefile}
  728. in the directory containing the sequences. You should then change the compiled
  729. model file name, by changing the value of the
  730. @code{model_file_name} variable to the name of your choice
  731. in the Makefile. It is set to @file{mymodel} in the template. For the
  732. predator-prey model, it could be set like
  733. @example
  734. model_file_name = predator
  735. @end example
  736. If you want the executable model file to be built in another directory, you could
  737. set
  738. @example
  739. model_file_name = some_dir/predator
  740. @end example
  741. The other items set in the default Makefile should be right.
  742. @end enumerate
  743. The preprocessing and the compilation are launched with
  744. @example
  745. make all
  746. @end example
  747. The mortran files are generated by the cmz directive preprocessor
  748. from files found in the package source directories. The mortran files
  749. end with @samp{.mtn} for the main files and @samp{.mti} for
  750. include files. They are output in the current directory.
  751. The mortran preprocessor then preprocess these mortran files and includes
  752. the sequences. The resulting fortran code is also in the current directory,
  753. in files with a @samp{.f} suffix.
  754. Some fortran files ending with @samp{.F} may also be
  755. created by the cmz directive preprocessor.
  756. The object files resulting from the compilation of all the
  757. fortran files (generated from mortran or directly from fortran files) are
  758. there too.
  759. In case you want to override the default sequences or a subroutine file
  760. you just have to create it in your working directory along with the
  761. @file{zinit.mti}. @c and @file{dimetaphi.mti}.
  762. For example you could want to
  763. create or modify a @file{zsteer.mti} file (@pxref{End of time step,,
  764. Executing code at the end of each time step}), a @file{zcmd_law.mti} file
  765. (@pxref{Control laws}), a @file{monitor.f} file
  766. (@pxref{Turning the model into a subroutine}) to take advantage of
  767. features presented later in this manual.
  768. More in-depth discussion of using make to run @Minik{} is covered in
  769. @ref{Advanced use of @Minik{} with make}.
  770. For example it is also possible to create files that are to be
  771. preprocessed by the cmz directive
  772. preprocessor and separate source files and generated files.
  773. This advanced use is more precisely covered in
  774. @ref{Programming with cmz directives}.
  775. @page
  776. @node Simulation and output
  777. @subsection Running a simulation and using the output
  778. @cindex running model
  779. Once compiled the model is ready to run, it only has to be executed. On
  780. standard output informations about the states, transfers, tangent linear
  781. system and other jacobian matrices are printed.
  782. For example the predator-prey model could be executed with:
  783. @example
  784. ./predator > result.lis
  785. @end example
  786. @cindex output file
  787. @vindex dEta(.)
  788. @cindex @file{res.data}
  789. @cindex @file{dres.data}
  790. @cindex @file{tr.data}
  791. @cindex @file{aspha.data}
  792. @cindex @file{Model.hlp}
  793. @c In case of a model entered symbolically
  794. @c (@pxref{Symbolic model description})
  795. The correspondance
  796. between the symbolic variables and the basic vectors and functions
  797. are printed at run time:
  798. @example
  799. ---------------- Informing on Phi definition -----------------
  800. Var-name, Function-name, index in ff vector
  801. ff_interact f_interact 1
  802. ----------------------------------------------------
  803. --------------- Informing on Eta definition ------------------
  804. Var-name, Function-name, index in eta vector
  805. eta_prey deta_prey 1
  806. eta_pred deta_pred 2
  807. @end example
  808. A summary of the model equations are in @file{Model.hlp} file. For
  809. the same example:
  810. @example
  811. ======================= set_Phi
  812. 1 ff_interact f_interact eta_pray*eta_pred
  813. ======================= set_Eta
  814. 1 eta_pray deta_pray apar*eta_pray-apar*ff_interact
  815. 2 eta_pred deta_pred -cpar*eta_pred+cpar*ff_interact
  816. @end example
  817. @c FIXME never talked about that. Certainly not here
  818. when other general functions are specified with @code{f_set}, it can appear
  819. also in the same help file when replaced by @code{fun_set}.
  820. As far as possible, all data printed in the listing are associated
  821. with a name related to a variable. Here is an extract:
  822. @example
  823. Gamma :-8.19100E-02-1.42151E-01 3.87150E-02
  824. eta_courant eta_T_czcx eta_T_sz
  825. ------------------------------------------------
  826. Omega : 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
  827. courant_L T_czcx Psi_Tczc Psi_Tsz
  828. ------------------------------------------------
  829. @end example
  830. for the two known vectors of the system, and:
  831. @example
  832. >ker : Matrice de couplage 4 4 4 4
  833. courant_L Raw(1,j=1,4): 1.000 -9.9010E-03 0.000 0.000
  834. T_czcx Raw(2,j=1,4): -2.7972E-02 1.000 0.000 9.9900E-04
  835. Psi_Tczcx Raw(3,j=1,4): 0.1605 9.7359E-02 1.000 -5.7321E-03
  836. Psi_Tsz Raw(4,j=1,4): 0.000 -0.1376 5.7225E-03 1.000
  837. Var-Name courant_L T_czcx Psi_Tczc Psi_Tsz
  838. ----------------------------------------------------------
  839. @end example
  840. where the @code{couplage} (coupling matrix) is given that corresponds
  841. to the matrix coupling the four transfer components after @math{\delta\eta}
  842. has been eliminated from system. It is computed in the subprogram
  843. @file{oker} (for kernel) which solves the system.
  844. Basic results are output in a set of @samp{.data} files.
  845. The first line (or two lines) describes the column with a @samp{#}
  846. character used to mark the lines as comments (for @command{gnuplot}
  847. for example).
  848. In the @samp{.data} files, the data are simply separated with spaces.
  849. Each data file has the @code{time} variable values as first column.
  850. @footnote{@file{dres.data} has another time related variable as second column:
  851. @cindex @file{dres.data}
  852. @vindex dt
  853. @code{dt}, the time step that can vary in the course of a simulation.}.
  854. Following columns give the values of @code{eta(.)} in @file{res.data},
  855. @code{dEta(.)} in @file{dres.data} -- the step by step variation of
  856. @code{eta(.)} -- and @code{ff(.)} in @file{tr.data}.
  857. Along the simulation the @acronym{TEF} Jacobian matrices are computed.
  858. A transfer variables elimination process also leads to the definition
  859. of the classical state advance matrix of the system
  860. (the corresponding array is @code{aspha(.,.)} in the code).
  861. This matrix is output in the file @file{aspha.data} that is used to
  862. post-run dynamics analyses. The matrix columns are written column wise on each
  863. record.
  864. @xref{Stability of fastest modes,,Stability analysis of fastest modes}.
  865. @xref{Generalized TLS,,Generalized
  866. tangent linear system analysis}. It is not used in the solving process.
  867. Other @samp{.data} files will be described later.
  868. @c FIXME already said
  869. @c At the beginning of a run, the help file @file{Model.hlp} is generated for
  870. @c global checkiing of the model. In the example, it is:
  871. @c @example
  872. @c ======================= set_Phi
  873. @c 1 ff_interact f_interact eta_pray*eta_pred
  874. @c ======================= set_Eta
  875. @c 1 eta_pray deta_pray apar*eta_pray-apar*ff_interact
  876. @c 2 eta_pred deta_pred -cpar*eta_pred+cpar*ff_interact
  877. @c @end example
  878. @node Graphics
  879. @subsection Doing graphics
  880. @cindex graphics
  881. @cindex graphics with @command{gnuplot}
  882. @cindex graphics with @command{PAW}
  883. @c The format of the @samp{.data} files are coherent with GNU graphics, that is
  884. @c the data are simply separated with spaces.
  885. Since the data are simply separated with spaces, and comment lines
  886. begin with @samp{#}, the
  887. files can be vizualised with many programs.
  888. With @command{gnuplot}, for example, to plot @code{eta(@var{n})},
  889. the @command{gnuplot} statement could be:
  890. @example
  891. plot "res.data" using 1:(@var{n}+1)
  892. @end example
  893. The similar one for @code{ff(@var{n})}:
  894. @example
  895. plot "tr.data" using 1:(@var{n}+1)
  896. @end example
  897. For people using @command{PAW}, the CERN graphical computer code,
  898. @Minik{} prepares
  899. kumacs that allow to read process the @samp{.data} files in the form of
  900. @emph{n-tuples} (see the @cite{PAW manual} for more information).
  901. In that cas, the flag @code{sel paw} has to be gievn in the @file{selsequ.kumac}.
  902. The generated n-tuples are ready to use only
  903. for vector dimension of at most 10 (including the variable @code{time}).
  904. These kumacs are overwritten each time the model is run. Usaually, gnuplot has
  905. to be preferred, but when using surfaces and histograms, PAW is better.
  906. The @file{gains.f} (and @file{go.xqt} is provided as an example in the
  907. @Minik{} files.
  908. @node Controlling the run
  909. @section Controlling the run
  910. @cindex controlling the run
  911. It is possible to add code that will be executed at the end of each time
  912. step. It is also possible to specify which time step leads to a printout on
  913. standard output. For maximal control, the code running te model may be
  914. turned into a subroutine to be called from another fortran (or C)
  915. program, this possibility is covered in @ref{Calling the model code}.
  916. @menu
  917. * End of time step::
  918. * Controlling the printout and data output::
  919. @end menu
  920. @node End of time step
  921. @subsection Executing code at the end of each time step
  922. @cindex @file{zsteer}
  923. @cindex @file{zsteer.inc}
  924. The code in the sequence @file{zsteer} is executed at the end of each time
  925. step. It is possible to change the time step length (variable @code{dt})
  926. verify that the non linearity are not too big, or perform
  927. discontinuous modifications of the states. One available variable @code{res}
  928. might be usefull for time step monitoring. At the end of the time step,
  929. as soon as @math{\varphi} has been computed, a numerical test is applied
  930. on a pseudo relative quadratic residual between
  931. @math{\varphi=f(\eta(t-dt)+d\varphi} (@code{ ffl}), where @math{d\varphi}
  932. is given by the system resolution in @code{ker},and
  933. @math{\varphi=f(\eta),\varphi)}, Fortran variable (@code{ff}):
  934. @verbatim
  935. ! ========================================================
  936. ! test linearite ffl - ff
  937. ! ========================================================
  938. if (istep.gt.1)
  939. < res=0.; <io=1,m; res = res +(ffl(io)-ff(io))**2/max(one,ff(io)*ff(io)); >;
  940. if (res .gt. TOL_FFL)
  941. < print*,'*** pb linearite : res > TOL_FFL a istep',istep,res,' > ',TOL_FFL;
  942. do io=1,m < z_pr: io,ff(io),ff(io)-ffl(io); >;
  943. >;
  944. >;
  945. @end verbatim
  946. This test hence applies only for non linearities in tranfer models. Nevertheless,
  947. @code{res} might be usefull to monitor the time step @code{dt} in @code{ZSTEER}
  948. and eventually go backward one step (@code{goto :ReDoStep:}).
  949. This can more appropriatly be coded in the (empty in default case)
  950. sequence @code{zstep}, inserted just before time-advancing
  951. states and @code{time} variables in @file{principal}.
  952. @vindex ffl(.)
  953. @cindex @code{ffl} (linearity test)
  954. @cindex linearity test
  955. It is also possible to fix the value of the criterium @code{TOL_FFL} in
  956. @file{zinit} different from its default value of @math{10^{-3}} --
  957. independent of the Fortran precision.
  958. Many other variables are available, including
  959. @vtable @code
  960. @item istep
  961. The step number;
  962. @item couplage(.)
  963. The @acronym{TEF} coupling matrix between transfers;
  964. @item H
  965. The Jacobian matrix corresponding with:
  966. @c \varphi(t) &= f(\eta(t),\varphi(t))\cr
  967. @c \frac{\partial g(\eta(t),\varphi(t))}{\partial \eta(t)}
  968. @tex
  969. $$\partial_{\eta} g(\eta(t),\varphi(t));
  970. $$
  971. @end tex
  972. @ifnottex
  973. g_1(eta,phi);
  974. @end ifnottex
  975. @item Bb
  976. The Jacobian matrix corresponding with:
  977. @tex
  978. $$\partial_{\varphi} g(\eta(t),\varphi(t));
  979. $$
  980. @end tex
  981. @ifnottex
  982. g_2(eta,phi);
  983. @end ifnottex
  984. @item Bt
  985. The Jacobian matrix corresponding with:
  986. @tex
  987. $$\partial_{\eta} f(\eta(t),\varphi(t));
  988. $$
  989. @end tex
  990. @ifnottex
  991. f_1(eta,phi);
  992. @end ifnottex
  993. @item D
  994. The Jacobian matrix corresponding with:
  995. @tex
  996. $$\partial_{\varphi} f(\eta(t),\varphi(t));
  997. $$
  998. @end tex
  999. @ifnottex
  1000. f_2(eta,phi);
  1001. @end ifnottex
  1002. @item aspha
  1003. The state advance matrix;
  1004. @item dneta
  1005. @itemx dphi
  1006. the variable increments;
  1007. @end vtable
  1008. One should be aware of that the linearity test concerns the preceding step.
  1009. We have yet no example of managing the time-step.
  1010. @page
  1011. @node Controlling the printout and data output
  1012. @subsection Controlling the printout and data output
  1013. @cindex printing
  1014. @vindex zprint
  1015. @vindex modzprint
  1016. The printout on standard output is performed if the variable @code{zprint}
  1017. of type @code{logical} is true. Therefore it is possible to control this
  1018. printout by setting @code{zprint} false or true. For example the following
  1019. code, in sequence @file{zsteer}, triggers printing for every
  1020. @code{modzprint} time step and the two following time steps:
  1021. @example
  1022. ZPRINT = mod(istep+1,modzprint).eq.0;
  1023. Zprint = zprint .or. mod(istep+1,modzprint).eq.1;
  1024. Zprint = zprint .or. mod(istep+1,modzprint).eq.2;
  1025. @end example
  1026. The data output to @file{.data} files described in @ref{Simulation and output,,
  1027. Running a simulation and using the output} is performed if the
  1028. @code{logical} variable @code{zout} is true. For example the following
  1029. code, in @file{zsteer}, triggers output to @file{.data} files every
  1030. @code{modzout} step.
  1031. @example
  1032. Zout = mod(istep,modzout).eq.0;
  1033. @end example
  1034. @node Advanced programming
  1035. @chapter Advanced @Minik{} programming
  1036. @menu
  1037. * Selecting features::
  1038. * Calling the model code::
  1039. * 1D gridded model::
  1040. * Double precision::
  1041. * Partial Derivatives::
  1042. * Rule of programming non continuous models::
  1043. * Parameters::
  1044. * Observations and data::
  1045. * Explicit model size::
  1046. * Programming with cmz directives::
  1047. @end menu
  1048. @node Selecting features
  1049. @section Overview of additional features setting
  1050. @cindex feature setting
  1051. @cindex select flag
  1052. @cindex logical flags
  1053. @cindex @file{selseq.kumac}
  1054. It is possible to enable some features by selecting which code should
  1055. be part of the principal program. Each of these optionnal features are
  1056. associated with a @dfn{select flag}.
  1057. For example
  1058. @c the optimisation with minuit is associated with the select
  1059. @c flag @samp{minuik},
  1060. double precision is used instead of simple precision with
  1061. the @samp{double} select flag,
  1062. the model is a subroutine with the select flag @samp{monitor},
  1063. the Kalman filter code is set with @samp{kalman} and the 1D gridded
  1064. model capabilities are associated with @samp{grid1d}.
  1065. @c Currently it is only possible
  1066. @c to select features in cmz.
  1067. To select a given feature the cmz statement @code{sel @var{select_flag}} should
  1068. be written down in the @file{selseq.kumac} found in the model directory.
  1069. With make either the corresponding variable should be set to 1 or it
  1070. should be added to the @code{SEL} make variable, depending on the feature.
  1071. Other features don't need different or additional code to be used.
  1072. Most of the features are enabled by setting specific logical variables to
  1073. @samp{.true.}. This is the case for
  1074. @code{zback} for the adjoint model, @code{zcommand} if the command is in a file
  1075. and @code{zlaw} if it is a function and @code{zkalman} for the Kalman filter.
  1076. These select and logical flags are described in the corresponding sections.
  1077. In cmz an alternative of writing select flags to @file{selseq.kumac} is to
  1078. drive the compilation with @code{smod @var{sel_flag}}. In that case the
  1079. @var{sel_flag} is selected and the files and executable goes to a directory
  1080. named @file{sel_flag}.
  1081. The select flags are taken into account during cmz directives preprocessing.
  1082. Therefore you have the possibility to use these flags to conditionnaly
  1083. include pieces of code. In most cases you don't need to include code conditionally
  1084. yourself though, but if you want to, this is covered in
  1085. @ref{Programming with cmz directives}.
  1086. @node Calling the model code
  1087. @section Calling the model code
  1088. When the model code is a subroutine, it can be called from another fortran
  1089. program unit (or another program), and the model will be
  1090. run each time the subroutine is called.
  1091. This technique could be used, for example to perform optimization
  1092. (@pxref{Adjoint model and optimisation,,Adjoint model and optimisation
  1093. with @Minik{}}), or to run the model with different parameters.
  1094. @menu
  1095. * Turning the model into a subroutine::
  1096. * The model subroutine::
  1097. @end menu
  1098. @node Turning the model into a subroutine
  1099. @subsection Turning the model into a subroutine
  1100. @c This feature is allready enabled with @command{make}, and to
  1101. @c override the default program a file called @file{monitor.f} has to be created
  1102. @c in the working directory. The default program simple calls the model
  1103. @c subroutine.
  1104. With cmz, one has to do a
  1105. @example
  1106. sel monitor
  1107. @end example
  1108. in the @file{selseq.kumac} file and create the KEEP that call the
  1109. model code. @xref{Selecting features, Overview of additional features
  1110. setting}.
  1111. With make @samp{monitor} should be added to the @code{SEL} variable in
  1112. the @file{Makefile}, for example:
  1113. @example
  1114. SEL = monitor
  1115. @end example
  1116. A file that call the principal subroutine should also be written, using
  1117. the prefered language of the user. The additional object files should
  1118. then be linked with the @Minik{} objects. To that aim they may be added
  1119. to the @code{miniker_user_objects} variable.
  1120. @node The model subroutine
  1121. @subsection Calling the model subroutine
  1122. The model subroutine is called @samp{principal} and is called with the
  1123. following arguments:
  1124. @deffn Subroutine principal (Cost, ncall, integer_flag, file_suffix, info, idxerror)
  1125. Where @var{Cost} is a real number, @code{real} or @code{double precision},
  1126. and is set by the @code{principal}
  1127. subroutine. It holds the value of the cost function if such function has been
  1128. defined (the use and setting of a cost function is covered later,
  1129. see @ref{Cost function coding and adjoint modeling}).
  1130. @var{ncall} is an integer which corresponds with the number of
  1131. call to @code{principal} done so far, it should be initialized to 0 and
  1132. its value should not be changed, as it is changed in the @code{principal}
  1133. subroutine.
  1134. @var{integer_flag} is an integer that can be set by the user to be accessed
  1135. in the @code{principal} subroutine. For example its value could be used to
  1136. set some flags in the @file{zinit} sequence.
  1137. @var{file_suffix} is a character string, that is suffixed to the output files
  1138. names instead of @samp{.data}. If the first character is the null character
  1139. @samp{char(0)}, the default suffix, @samp{.data} is appended.
  1140. @var{info} and @var{idxerror} are integer used for error reporting.
  1141. @var{idxerror} value is 0 if there was no error. It is negative for
  1142. an alert, positive for a very serious error. The precise value determines
  1143. where the error occured.
  1144. @var{info} is an integer holding more precise information about the
  1145. error. It is usually the information value from lapack.
  1146. The precise meaning of these error codes is in @ref{tab:error_codes}.
  1147. @end deffn
  1148. @float table, tab:error_codes
  1149. @multitable {kalman analysis state matrix advance in phase space, @math{(I-D)} inversion} {inversion} {@code{idxerror}}
  1150. @headitem Source of error or warning @tab @code{info} @tab @code{idxerror}
  1151. @c @item @code{} @tab @file{.data} @tab @tab
  1152. @item state matrix inversion in ker @tab inversion @tab 1
  1153. @item time advance system resolution in ker @tab system @tab 2
  1154. @item transfer propagator, @math{(I-D)} inversion @tab inversion @tab 3
  1155. @item kalman analysis state matrix advance in phase space, @math{(I-D)} inversion @tab inversion @tab 21
  1156. @item kalman analysis variance covariance matrix non positive @tab Choleski @tab 22
  1157. @item kalman analysis error matrix inversion @tab inversion @tab 23
  1158. @item kalman error matrix advance @tab system @tab 24
  1159. @item transfers determination linearity problem for transfers @tab @tab -1
  1160. @item transerts determination Newton D_loop does not converge @tab @tab -2
  1161. @end multitable
  1162. @caption{Meaning of error codes returned by principal.}
  1163. @end float
  1164. In general more information than the provided arguments has to be passed
  1165. to the @code{principal} subroutine, in that case a @code{common} block,
  1166. to be written in the @file{zinit} sequence can be used.
  1167. @page
  1168. @node 1D gridded model
  1169. @section Describing 1D gridded model
  1170. Specific macros have been built that allow generic description of 1D gridded models.
  1171. Because of the necessity of defining left and right limiting conditions, the models
  1172. are partitionned in three groups for cell and transfer components. In the following example,
  1173. a chain of masselottes linked by springs and dumps is bounded to a wall on the left,
  1174. and open at right. The @acronym{TEF} formulation of the problem is written in the phase space (position-shift, velocity)
  1175. for node @math{k}, with bounding conditions:
  1176. @ifset latex
  1177. @tex
  1178. \begin{equation}
  1179. \left\{
  1180. \begin{array}{cc}
  1181. \partial_t \eta _{k} ^{pos} = \eta _{k} ^{vel}\qquad& \\
  1182. \partial_t \eta _{k} ^{vel} = ( \varphi_k ^{spr} -\varphi _{k+1} ^{spr}&+\,\,\varphi _{k} ^{dmp}-\varphi _{k+1} ^{dmp})\,/m_k \nonumber\\
  1183. \end{array}
  1184. \right.
  1185. \end{equation}
  1186. \begin{equation}
  1187. \left\{
  1188. \begin{array}{cc}
  1189. \varphi_k ^{spr} = -k_k (\eta _{k} ^{pos}- \eta _{k-1} ^{pos})\\
  1190. \varphi_k ^{spr} = -d_k (\eta _{k} ^{vel}- \eta _{k-1} ^{vel})
  1191. \end{array}
  1192. \right.
  1193. \label{mass}
  1194. \end{equation}
  1195. \begin{equation}
  1196. \left\{
  1197. \begin{array}{cc}
  1198. \eta ^{pos}_{0} =& 0\\
  1199. \eta ^{vel}_{0} =& 0\\
  1200. \varphi ^{spr}_{N+1} =& 0\\
  1201. \varphi ^{dmp}_{N+1} =& 0
  1202. \end{array}
  1203. \right.
  1204. \end{equation}
  1205. @end tex
  1206. @end ifset
  1207. @ifclear latex
  1208. @tex
  1209. $$\left\{\eqalign{\partial_t \eta _{k} ^{pos} &= \eta _{k} ^{vel} \cr
  1210. \partial_t \eta _{k} ^{vel} &= ( \varphi_k ^{spr} -\varphi _{k+1} ^{spr} + \varphi _{k} ^{dmp}-\varphi _{k+1} ^{dmp})\,/m_k \cr}\right.$$
  1211. $$\left\{\eqalign{
  1212. \varphi_k ^{spr} &= -k_k (\eta _{k} ^{pos}- \eta _{k-1} ^{pos})\cr
  1213. \varphi_k ^{spr} &= -d_k (\eta _{k} ^{vel}- \eta _{k-1} ^{vel})
  1214. \cr}\right.$$
  1215. $$\left\{\eqalign{\eta ^{pos}_{0} &= 0\cr
  1216. \eta ^{vel}_{0} &= 0\cr
  1217. \varphi ^{spr}_{N+1} &= 0\cr
  1218. \varphi ^{dmp}_{N+1} &= 0\cr}\right.$$
  1219. @end tex
  1220. @ifnottex
  1221. States:@*
  1222. @noindent @math{d position(t,k)/d t = velocity(t,k)@*
  1223. d velocity (t,k)/d t =(spring(t,k) - spring(t,k+1)+ dump(t,k)- dump(t,k+1))/m_k}
  1224. Transfers:@*
  1225. @noindent @math{spring(t,k) = -k_k (position(t,k)- position(t,k-1))@*
  1226. dump(k,t) &= -d_k (velocity(t,k)- velocity(t,k-1))}
  1227. Bounding conditions:@*
  1228. @noindent @math{position(t,0) = 0@*
  1229. velocity(t,0) = 0@*
  1230. spring(t,N+1) = 0@*
  1231. dump(t,N+1) =0}
  1232. @end ifnottex
  1233. @end ifclear
  1234. @cindex down node
  1235. @cindex up node
  1236. where @math{m_k} is the mass of node @math{k}, @math{r_k} and @math{d_k}
  1237. the rigidity of springs and dumping coefficients. There are @math{N} nodes
  1238. in the grid, from 1 to @math{N}, and two nodes outside of the grid,
  1239. a limiting node 0, and a limiting node @math{N+1}. The limiting node
  1240. corresponding with node 0 is called the @dfn{down} node, while the limiting
  1241. node corresponding with node @math{N+1} is called the @dfn{up} node.
  1242. Other models not part of the 1D grid may be added if any.
  1243. To enable 1D gridded models, one should set the select flag @samp{grid1d}.
  1244. In cmz it is achieved setting the select flag in
  1245. @file{selseq.kumac}, like
  1246. @example
  1247. sel grid1d
  1248. @end example
  1249. With make, the @code{SEL} variable should contain @code{grid1d}. For example
  1250. to select @code{grid1d} and @code{monitor}, it could be
  1251. @example
  1252. SEL = grid1d,monitor
  1253. @end example
  1254. @menu
  1255. * 1D gridded Model size::
  1256. * 1D gridded model code::
  1257. @end menu
  1258. @node 1D gridded Model size
  1259. @subsection Setting dimensions for 1D gridded model
  1260. @c FIXME grid1d sans dimetaphi?
  1261. In that case the number of nodes, the number of states and tranferts
  1262. per node, and the number of limiting transfers and states are required.
  1263. These dimensions has to be entered in the
  1264. @file{DimEtaPhi} sequence. The parameters for cells are
  1265. @vtable @code
  1266. @item n_node
  1267. Number of cell nodes in the 1D grid.
  1268. @item n_dwn
  1269. Number of limiting cells with index -1, @i{i.e.} number of cells in the
  1270. limiting down node.
  1271. @item n_up
  1272. Number of limiting cells with index +1, @i{i.e.} number of cells in the
  1273. limiting up node.
  1274. @item n_mult
  1275. Number of cells in each node (multiplicity).
  1276. @end vtable
  1277. @vindex m_node
  1278. @vindex m_dwn
  1279. @vindex m_up
  1280. @vindex m_mult
  1281. The parameters for transfers, are similarly
  1282. @code{m_node}, @code{m_dwn}, @code{m_up}, @code{m_mult}.
  1283. The layout of their declaration should be respected as
  1284. the precompiler matches the line. Also this procedure is tedious, it
  1285. should be selected for debuging processes (use the flag @code{sel dimetaphi}
  1286. in ``selsequ.kumac''. Otherwise, the dimensioning sequence will be automaticaly
  1287. generated, which is smart but can lead to diffculty in interpreting syntax errors.
  1288. Once a model is correctly entred, turn off the sel flag and further modifications
  1289. will automatically generate the proper dimensions. The correctness of dimensionning
  1290. should nevertheless always be checked in @code{principal.f}, where you can also
  1291. check that null valued parameters as @code{lp, mobs, nxp} will suppress parts
  1292. of the code - this is signaled as Fortran comment cards.
  1293. In our example, there are three grids of cell and
  1294. transfer variables (@code{n_node=m_node=3}).
  1295. There are two cells and two transfers in each node
  1296. (@code{n_mult=2} and @code{m_mult=2}). There is no limiting condition
  1297. for the states in the down node therefore @code{n_up=0}.
  1298. There is no transfer for the first limiting node, and
  1299. therefore @code{m_dwn=0}.
  1300. There are two states in the limiting node 0, the down node,
  1301. @code{n_dwn=2}, and two transfers in the limiting last node the node up,
  1302. and @code{m_up=2}:
  1303. @example
  1304. ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  1305. ! nodes parameters, and Limiting Conditions (Low and High)
  1306. ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  1307. parameter (n_node=3,n_dwn=2,n_up=0,n_mult=2);
  1308. parameter (m_node=3,m_dwn=0,m_up=2,m_mult=2);
  1309. ! ________________________________________________________
  1310. @end example
  1311. @ignore
  1312. @c FIXME enleve par al1
  1313. The dimension of the parameter arrays should also be declared in the
  1314. @file{dimetaphi} sequence. Here we have 3 parameters, for
  1315. @math{m_k}, @math{r_k} and @math{d_k}:
  1316. @example
  1317. dimension rk(n_node),rd(n_node),rmassm1(n_node);
  1318. @end example
  1319. @end ignore
  1320. @node 1D gridded model code
  1321. @subsection 1D gridded Model coding
  1322. The model code and parameters go in the @file{zinit} sequence.
  1323. @subsubheading Parameters
  1324. A value for the @Minik{} parameters and the model parameters should be given in
  1325. @file{zinit}, in our example we have
  1326. @example
  1327. !%%%%%%%%%%%%%%%%%%%%%%
  1328. ! Parameters
  1329. !%%%%%%%%%%%%%%%%%%%%%%
  1330. real rk(n_node),rd(n_node),rmassm1(n_node);
  1331. data rk/n_node*1./;
  1332. data rd/n_node*0.1/;
  1333. data rmassm1/n_node*1./;
  1334. dt=.01;
  1335. nstep=5 000;
  1336. modzprint = 1000;
  1337. time=0.;
  1338. @end example
  1339. @subsubheading Limiting conditions
  1340. @cindex limiting conditions
  1341. @c The limiting states and transfer variables and the corresponding equations are
  1342. @c declared using
  1343. @c the symbolic model description
  1344. @c (@pxref{Symbolic model description}).
  1345. There are four mortran blocks for @code{node} and @code{up} and @code{down}, both
  1346. for states and transfers:
  1347. @findex set_dwn_eta
  1348. @findex set_dwn_phi
  1349. @findex set_up_eta
  1350. @findex set_up_phi
  1351. @table @code
  1352. @item set_dwn_eta
  1353. down node cells
  1354. @item set_up_eta
  1355. up node cells
  1356. @item set_dwn_phi
  1357. down node transfers
  1358. @item set_up_phi
  1359. up node transfers
  1360. @end table
  1361. The following scheme illustrates the example:
  1362. @smallexample
  1363. !%%%%%%%%%%%%%%%%%%%%%%%%%%================================================
  1364. ! Maillage convention inode
  1365. !%%%%%%%%%%%%%%%%%%%%%%%%%% Open ended
  1366. !(2 Down Phi Eta (n_node)
  1367. ! Eta) \| .-----. .-----. .-----. /
  1368. ! wall \|-\/\/\-| |-\/\/\-| | . . . -| |-\/\/\- |dummy
  1369. ! pos \|--***--| 1 |--***--| 2 | . . . -| n |--***-- |Phis
  1370. ! speed \| 1 |_____| 2 |_____| n |_____| n+1 \(2 Up Phi)
  1371. !
  1372. @end smallexample
  1373. Two states are associated with the down node, they correspond to the position
  1374. and speed of the wall. As the wall don't move these states are initialized to
  1375. be 0, and the cells are stationnary cells, therefore these values remain 0.
  1376. @example
  1377. ! Down cells (wall)
  1378. ! -----------------
  1379. eta_pos_wall = 0; eta_speed_wall = 0.;
  1380. set_dwn_eta
  1381. < var: eta_pos_wall, fun: deta_pos_wall = 0.;
  1382. var: eta_speed_wall, fun: deta_speed_wall= 0.;
  1383. >;
  1384. @end example
  1385. There are 2 limiting transfers in the up node. They correspond with an open
  1386. end and are therefore set to 0.
  1387. @example
  1388. ! limiting Transfers : dummy ones
  1389. ! -------------------------------
  1390. set_Up_Phi
  1391. < var:ff_dummy_1, fun: f_dummy_1=0.;
  1392. var:ff_dummy_2, fun: f_dummy_2=0.;
  1393. >;
  1394. @end example
  1395. @subsubheading Starting points
  1396. The cell node state values are initialized. They are in an array
  1397. indexed by the @code{inode} variable. In the example the variable
  1398. corresponding with position is @code{eta_move} and the variable corresponding
  1399. with speed is @code{eta_speed}. Their initial values are set with the
  1400. following mortran code
  1401. @example
  1402. !---------------
  1403. ! Initialisation
  1404. !---------------
  1405. ;
  1406. do inode=1,n_node <eta_move(inode)=0.01; eta_speed(inode)=0.0;>;
  1407. @end example
  1408. If any transfer needs to be given a first-guess value, this is also done
  1409. using @code{inode} as the node index.
  1410. @subsubheading Grid node equations
  1411. @findex set_node_Phi
  1412. @findex set_node_eta
  1413. @cindex equations, grid
  1414. Each node is associated with an index @code{inode}. It allows to refer to the
  1415. preceding node, with @code{inode-1} and the following node @code{inode+1}.
  1416. The node states are declared in @code{set_node_Eta} block and the transfers are
  1417. in @code{set_node_Phi} blocks.
  1418. In the example, the cells are declared with
  1419. @example
  1420. ! node cells
  1421. ! ----------
  1422. ;
  1423. set_node_Eta
  1424. < var: eta_move(inode), fun: deta_move(inode) = eta_speed(inode);
  1425. var: eta_speed(inode),
  1426. fun: deta_speed(inode) = rmassm1(inode)
  1427. *( - ff_spring(inode+1) + ff_spring(inode)
  1428. - ff_dump(inode+1) + ff_dump(inode)
  1429. );
  1430. >;
  1431. @end example
  1432. Note that the @code{inode} is dummy in the @code{var:} definition and can as
  1433. well be written as: @code{var: eta_move(.)}.
  1434. The transfers are (@code{ff_spring} corresponds with springs and
  1435. @code{ff_dump} with dumps):
  1436. @example
  1437. !%%%%%%%%%%%%%%%%%%%%%%
  1438. ! Transfer definition
  1439. !%%%%%%%%%%%%%%%%%%%%%%
  1440. ! node transfers
  1441. ! --------------
  1442. ! convention de signe spring : comprime:= +
  1443. set_node_Phi
  1444. < var: ff_spring(.),
  1445. fun:
  1446. f_spring(inode)= -rk(inode)*(eta_move(inode) - eta_move(inode-1));
  1447. var: ff_dump(.),
  1448. fun:
  1449. f_dump(inode) = -rd(inode)*(eta_speed(inode) - eta_speed(inode-1));
  1450. >;
  1451. @end example
  1452. The limiting states and transfers are associated with the states or transfers
  1453. with index @code{inode+1} or @code{inode-1} appearing in node cell and
  1454. transfer equations (@code{inode-1} for down limiting conditions and
  1455. @code{inode+1} for up limiting conditions) in their order of appearance.
  1456. In our example, in the @code{eta_speed} state node equation
  1457. @code{ff_spring(inode+1)} appears before @code{ff_dump(inode+1)} and is
  1458. therefore associated with @code{ff_dummy_1} while @code{ff_dump(inode+1)} is
  1459. associated with the @code{ff_dummy_2} limiting transfer, as @code{ff_dummy_1}
  1460. appears before @code{ff_dummy_2} in the limiting up transfers definitions.
  1461. Verification of the grid index coherence should be eased with the following
  1462. help printed in the listing header:
  1463. @example
  1464. --------------- Informing on Dwn Eta definition ---------------
  1465. Var-name, Function-name, index in eta vector
  1466. eta_pos_wall deta_pos_wall 1 [
  1467. eta_speed_wall deta_speed_wall 2 [
  1468. -------------- Informing on Eta Nodes definition --------------
  1469. Var-name, Function, k2index of (inode: 0 [ 1,...n_node ] n_node+1)
  1470. eta_move deta_move 1 [ 3 ... 7 ] 9
  1471. eta_speed deta_speed 2 [ 4 ... 8 ] 10
  1472. ---------------- Informing on Up Phi definition -------------
  1473. Var-name, Function-name, index in ff vector
  1474. ff_dummy_1 f_dummy_1 ] 7
  1475. ff_dummy_2 f_dummy_2 ] 8
  1476. ff_move_sum f_move_sum ] 9
  1477. ff_speed_sum f_speed_sum ] 10
  1478. ----------------------------------------------------
  1479. -------------- Informing on Phi Nodes definition ---------------
  1480. Var-name, Function, k2index of (inode: 0 [ 1,...m_node ] m_node+1)
  1481. ff_spring f_spring -1 [ 1 ... 5 ] 7
  1482. ff_dump f_dump 0 [ 2 ... 6 ] 8
  1483. ----------------------------------------------------
  1484. @end example
  1485. All variable names and functions are free but has to be
  1486. different.
  1487. Any particular node-attached variable @math{k} is referred to as: @samp{(inode:k)},
  1488. where @math{k} has to be a Fortran expression allowed in arguments. The symbol
  1489. @samp{inode} is
  1490. reserved. As usual other Fortran instructions can be written within the
  1491. Mortran block @samp{< >} of each @code{set_} block.
  1492. @node Double precision
  1493. @section Double precision
  1494. The default for real variables is the @code{real} Fortran type. It is possible to
  1495. use double precision instead. In that case all the occurences of @samp{real@ }
  1496. in mortran code is substituted with @samp{double precision@ } at
  1497. precompilation stage,
  1498. and the Lapack subroutine names are replaced by the double precision names.
  1499. Eventual users'declaration of @code{complex@ } Fortran variables is also
  1500. changed to @code{double complex@ }.
  1501. This feature is turned on by @code{sel double} in @file{selseq.kumac} with cmz
  1502. and @code{double = 1} in the @file{Makefile} with make.
  1503. In order for the model to run as well in double as in simple precision,
  1504. some care should be taken to use the generic intrinsic functions, like
  1505. @code{sin} and not @code{dsin}. No numerical constant should be passed directly
  1506. to subroutines or functions, but instead a variable with the right type should
  1507. be used to hold the constant value, taking advantage of the implicit casts
  1508. to the variable type.
  1509. @node Partial Derivatives
  1510. @section Partial Derivatives
  1511. The partial derivative rules are included in a @code{Mortran} macro series
  1512. in @file{Derive_mac} of @Minik{} files. When using an anusual function,
  1513. one should verify that the corersponding rules are in that file.
  1514. It is easy to understand and add new rules in analogy with the already existing ones.
  1515. For instance, suppose one wants to use the intrinsic Fortran function @code{ abs()}.
  1516. Its derivatives uses the other function @code{sign()} this way:
  1517. @example
  1518. &'(ABS(#))(/#)' = '((#1)(/#2)*SIGN(1.,#1))'
  1519. @end example
  1520. In such cases when one is adding a new rule, it is important to use the generic function names
  1521. only (i.e. @code{sin} not @code{dsin}), because when compilating @Minik{} in the double precision
  1522. version, or complex version, the generic names will correctly handle the different variable
  1523. types - which is not the case when coding with specific function names.
  1524. @menu
  1525. * Derivating a power function::
  1526. @end menu
  1527. @node Derivating a power function
  1528. @subsection Derivating a power function
  1529. Partial derivative of a function in exponent is not secure in its Fortran form
  1530. @code{g(x,y)**(f(y))}. It should be replaced by @code{power(g,f)} of
  1531. the @Minik{} @file{mathlib},
  1532. or by the explicit form @code{exp(f(y)*log(g(x,y)))}.
  1533. Its derivative will have the following form:
  1534. @ifset latex @c ***PBAl1
  1535. @tex
  1536. \begin{equation}
  1537. \partial_x f^g=g f^{g-1}\partial_x f + f^g \log f\partial_x g = f^{g-1}(g\partial_x f + f\partial_x g)
  1538. \end{equation}
  1539. @end tex
  1540. @end ifset
  1541. @ifclear latex
  1542. @tex
  1543. $$\eqalign{\partial_x f^g &= g f^{g-1}\partial_x f + f^g \log f\partial_x g\cr
  1544. &= f^{g-1}(g\partial_x f + f\partial_x g)\cr}$$
  1545. @end tex
  1546. @end ifclear
  1547. and is in the macros list already defined in: @file{DERIVE_MAC}.
  1548. @node Rule of programming non continuous models
  1549. @section Rule of programming non continuous models
  1550. Some models may originally be non continuous, as the ones using a Fortran instruction @code{IF}.
  1551. Some may use implicitly a step function on a variable. In such cases, the model has to be
  1552. set in a derivable form, and use a ``smooth step'' instead.
  1553. One should be aware of that this apparently mathematical treatment currently
  1554. indeed leads to a physical question about the macroscopic form of a physical law.
  1555. At a macroscipic level, a step function is usually a nonsense.
  1556. @cindex Heaviside function
  1557. Taking
  1558. the example of phase-change, a fluid volume does not change phase at once, and a ``smooth
  1559. change of state'' is a correct macroscopic model.
  1560. @Minik{} provides with the smooth step function
  1561. @emph{Heavyside}@footnote{This naming is a joke
  1562. for ``Inert'' Heaviside function.} in the @Minik{} @file{mathlib}:
  1563. @example
  1564. Delta = -1."K";
  1565. A_Ice = heavyside("in:" (T_K-Tf), Delta, "out:" dAIce_dT);
  1566. @end example
  1567. in this example, @code{Tf} is the ice fusion-temperature, @code{A_ice}
  1568. gives the ice-fraction
  1569. of the mesh-volume of water at temperature @code{T_k}.
  1570. The smooth-step function is a quasi
  1571. hyperbolic tangent function of @math{x/\Delta},
  1572. normalised from 0 to 1, with a maximum slope
  1573. of 2.5, see figure @ref{heavy}.
  1574. @float Figure, heavy
  1575. @image{heavyside}
  1576. @caption{Heaviside function and derivative}
  1577. @end float
  1578. @ignore
  1579. @tex PBAl1
  1580. \begin{figure}[h]
  1581. \psfig{figure=heavyside.ps,%
  1582. bbllx=60pt,bblly=180pt,bburx=526pt,bbury=650pt,width=10cm,clip=}
  1583. \caption{La fonction ``domptée'' de Heaviside et sa dérivée pour une variable adimensionnée}
  1584. \label{heavy}
  1585. \end{figure}
  1586. @end tex
  1587. @end ignore
  1588. For @code{Mortran} to be able to symbolicaly compute the partial derivarives, the rule
  1589. is in the table of macros as:
  1590. @example
  1591. &'(HEAVYSIDE(#,#,#))(/#)' = '((#1)(/#4)*HEAVYDELTA(#1,#2,#3))'
  1592. @end example
  1593. which uses the Foratn entry point @code{HeavyDelta} in the Fortrsan function @code{heavyside}.
  1594. Another type of problem arises when coding a
  1595. @code{var=min(f(x),g(x))} Fortran instruction.
  1596. In such a case one does not want a derivative and one will code:
  1597. @example
  1598. var = HeavySide(f(x)-g(x),Delta,dum)*g(x) + (1.-HeavySide(f(x)-g(x),Delta,dum)*f(x);
  1599. @end example
  1600. or equivalently:
  1601. @example
  1602. var = HeavySide(f(x)-g(x),Delta,dum)*g(x) + HeavySide(g(x)-f(x),-Delta,dum)*f(x);
  1603. @end example
  1604. @strong{Warning}: the value of the argument @var{Delta} is important because
  1605. it will fix the maximum
  1606. slope of the function that will appear as a coefficient in the
  1607. Jacbian matrices.
  1608. @node Parameters
  1609. @section Parameters
  1610. It is possible to specify some Fortran variables as specific model parameters.
  1611. Model parameters
  1612. may be used in sensitivity studies (@pxref{Sensitivity to a parameter})
  1613. and in the adjoint model (@pxref{Sensitivity of cost function to parameters}).
  1614. Nothing special is done with parameters with Kalman filtering.
  1615. @findex Free_parameter
  1616. The parameters are fortran variables that should be initialized somewhere
  1617. in @file{zinit}. For a variable to be considered as a parameter, it should
  1618. be passed as an
  1619. argument to the @code{Free_parameters} macro. For example if
  1620. @code{apar} and @code{cpar} (from the predator example) are to be considered
  1621. as parameters, @code{Free_parameters} should be called with:
  1622. @example
  1623. Free_parameter: apar, cpar;
  1624. @end example
  1625. @c Forward sensitivities are explained later (@pxref{Sensitivity to a parameter}),
  1626. @c the syntax only is described here.
  1627. When used with grid1d models (@pxref{1D gridded model,,
  1628. Describing 1D gridded model}) the @code{inode} number may appear in
  1629. parenthesis:
  1630. @example
  1631. Free_parameter: rd(1), rk(2);
  1632. @end example
  1633. @node Observations and data
  1634. @section Observations and data
  1635. Some support for observations and interactions with data is available.
  1636. The observations are functions of the model variables. They don't have
  1637. any action on the model result, but they may (in theory) be observed
  1638. and measured. The natural use of these observations is to be compared
  1639. with data that correspond with the values from real measurements.
  1640. They are used in the Kalman filter (@pxref{Kalman filter}).
  1641. The (model) observation vector is noted @math{\omega}
  1642. @c FIXME is seems untrue?
  1643. @c in this section ($\mu$ elsewhere,
  1644. and the observation function is noted @math{h}:
  1645. @tex
  1646. $$
  1647. \omega = h ( \eta , \varphi)
  1648. $$
  1649. @end tex
  1650. @ifnottex
  1651. @noindent @math{omega(t) = h(eta(t), phi(t))}
  1652. @end ifnottex
  1653. @menu
  1654. * Observations::
  1655. * Data::
  1656. @end menu
  1657. @node Observations
  1658. @subsection Observations
  1659. @vindex mobs
  1660. The observation functions are set in a @code{set_probe} block in
  1661. the @file{zinit} sequence.
  1662. @cindex observation function
  1663. @c FIXME doesn't exist anymore
  1664. @c @defmac eqn: Obs_tef(@var{i}) = f(eta(.),ff(.))
  1665. @c This macro defines the observation equation as usual in a @code{set_block<}.
  1666. @c @code{f} is a fortran
  1667. @c expression which may be function of cell state variables,
  1668. @c @samp{eta(1)}@dots{}@samp{eta(np)} and transfers
  1669. @c @samp{ff(1)}@dots{}@samp{ff(mp)}, or of course their symbolic names.
  1670. @c @end defmac
  1671. For example suppose that, in the predator-prey model, we only
  1672. have access to the total population of preys and predators, we would have:
  1673. @example
  1674. set_probe
  1675. < eqn: pop = eta_pred + eta_pray;
  1676. >;
  1677. @end example
  1678. @c it is always turned on, now
  1679. @c The corresponding code is used with @code{sel obs} in @file{selseq.kumac}
  1680. @c with cmz and @code{obs = 1} in @file{Makefile} with make. And the feature
  1681. @c is turned on and off at run time with the logical flag @code{zobs} corresponding
  1682. @c to an available data from measurement
  1683. @c @vindex etaobs(.)
  1684. @cindex @file{obs.data}
  1685. The number of observations is put in the integer variable @code{mobs}.
  1686. The observation vector corresponds with the part of the @code{ff(.)}
  1687. array situated past the regular transferts, @code{ff(mp+.)}, and is output
  1688. in the file @file{obs.data}.
  1689. @c @vindex obetad(.,.)
  1690. @c @vindex obephid(.,.)
  1691. @c @vindex obspha(.,.)
  1692. @node Data
  1693. @subsection Data
  1694. @vindex zgetobs
  1695. @vindex vobs(.)
  1696. @cindex @file{data.data}
  1697. Currently this code is only used if the Kalman code is activated. This
  1698. may be changed in the future.
  1699. The convention for data is that whenever some data are available, the
  1700. logical variable @code{zgetobs} should be set to @samp{.true.}. And the
  1701. @code{vobs(.)} vector should be filled with the data values. This
  1702. vector has the same dimension than the observation
  1703. vector and each coordinate is meant to correspond with one
  1704. coordinate of the observation vector.
  1705. This feature is turned on by setting the logical variable @code{zdata}
  1706. to @samp{.true.}, and the @code{zgetobs} flag is typically set in the
  1707. @file{zsteer} sequence (@pxref{End of time step,,Executing code at
  1708. the end of each time step}).
  1709. Every instant data are available (@code{zgetobs} is true) the observations
  1710. are written to the file @file{data.data}. With the Kalman filter more
  1711. informations are output to the @file{data.data} file,
  1712. see @ref{Kalman filter results}.
  1713. @node Explicit model size
  1714. @section Entering model size explicitely
  1715. It is possible to enter the model dimensions explicitely, instead of
  1716. generating them automatically, as it was done previously.
  1717. This feature is turned on by @code{sel dimetaphi}
  1718. in @file{selseq.kumac} with cmz
  1719. and @code{dimetaphi} added to the @code{SEL} variable in
  1720. the @file{Makefile} with make.
  1721. @menu
  1722. * Size sequence::
  1723. * Model with explicit size::
  1724. @end menu
  1725. @node Size sequence
  1726. @subsection The explicit size sequence
  1727. @cindex dimetaphi
  1728. @cindex model size
  1729. @vindex np
  1730. @vindex mp
  1731. @vindex maxstep
  1732. @cindex @file{dimetaphi}
  1733. The dimension of the model is entered in the sequence @file{dimetaphi},
  1734. using the fortran @code{parameter np} for @code{eta(.)} and
  1735. @code{mp} for @code{ff(.)}.
  1736. For the Lotka-Volterra model, we have two cell components and only one transfer.
  1737. @example
  1738. parameter (np=2,mp=1);
  1739. @end example
  1740. You should not change the layout of the parameter statement as the
  1741. mortran preprocessor matches the line.
  1742. You also have to provide other parameters even if you don't have any
  1743. use for them. If you don't it will trigger fortran errors.
  1744. It includes the @code{maxstep} parameter that can have any value but 0,
  1745. @code{lp} and @code{mobs} that should be 0 in the example, and @code{nxp},
  1746. @code{nyp} and @code{nzp} that should also be 0.
  1747. The layout is the following:
  1748. @example
  1749. parameter (np=2,mp=1);
  1750. parameter (mobs=0);
  1751. parameter (nxp=0,nyp=0,nzp=0);
  1752. parameter (lp=0);
  1753. parameter (maxstep=1);
  1754. @end example
  1755. If there are observations, (@pxref{Observations}), the
  1756. size of the observation vector is set in the @file{dimetaphi} sequence
  1757. by the @code{mobs} parameter. For example if there is one observation:
  1758. @example
  1759. parameter (mobs=1);
  1760. @end example
  1761. To specify parameters (@pxref{Parameters}), the number of such parameters
  1762. has to be declared in @file{dimetaphi} with the parameter @code{lp}.
  1763. Then, if there are two parameters, they are first declared with
  1764. @example
  1765. parameter (lp=2);
  1766. @end example
  1767. @node Model with explicit size
  1768. @subsection Entering the model equations, with explicit sizes
  1769. @cindex model equations
  1770. @findex Phi_tef(.)
  1771. @findex deta_tef(.)
  1772. @vindex eta(.), explicit sizes
  1773. @vindex ff(.), explicit sizes
  1774. When sizes are explicit, another possibility exists for entering
  1775. the model equations. The use of symbolic names, as described in
  1776. @ref{Model equations} is still possible, and it also becomes possible to
  1777. set directly the equations associated with the @code{eta(.)}
  1778. and @code{ff(.)} vectors.
  1779. In case the symbolic names are not used,
  1780. the model equations for cells and transfers are entered using a mortran macro,
  1781. @code{f_set}@footnote{@code{fun_set}, or equivalently @code{f_set}, is a
  1782. general mortran macro associating a symbol with a fortran expression.
  1783. Here, it is the name of the symbol (@code{eta}) that has a particular meaning
  1784. for the building of the model.}, setting the @code{eta(.)} evolution with
  1785. @code{deta_tef(.)}
  1786. and the transfer definitions @code{ff(.)} with @code{Phi_tef(.)}.
  1787. @defmac f_set Phi_tef(@var{i}) = f(eta(.),ff(.))
  1788. This macro defines the transfer @var{i} static equation.
  1789. @code{f} is a fortran
  1790. expression which may be function of cell state variables,
  1791. @samp{eta(1)}@dots{}@samp{eta(np)} and transfers
  1792. @samp{ff(1)}@dots{}@samp{ff(mp)}.
  1793. @end defmac
  1794. In the case of the predator-prey model, the transfer definition for
  1795. @math{\varphi_{meet}} is:
  1796. @example
  1797. f_set Phi_tef(1) = eta(1)*eta(2);
  1798. @end example
  1799. @defmac f_set deta_tef(@var{i}) = g(eta(@var{i}),ff(.))
  1800. This macro defines the cell state component @var{i} time evolution model.
  1801. @code{g} is a expression which may be function of cell state variables,
  1802. @samp{eta(1)}@dots{}@samp{eta(np)} and transfers
  1803. @samp{ff(1)}@dots{}@samp{ff(mp)}.
  1804. @end defmac
  1805. The two cell equations of the predator-prey model are, with index 1 for the
  1806. prey (@math{\eta_{prey}}) and index 2 for the predator (@math{\eta_{pred}}):
  1807. @example
  1808. f_set deta_tef(1) = apar*eta(1)-apar*ff(1);
  1809. f_set deta_tef(2) = - cpar*eta(2) + cpar*ff(1);
  1810. @end example
  1811. The whole model is:
  1812. @example
  1813. !%%%%%%%%%%%%%%%%%%%%%%
  1814. ! Transfer definition
  1815. !%%%%%%%%%%%%%%%%%%%%%%
  1816. ! rencontres (meeting)
  1817. f_set Phi_tef(1) = eta(1)*eta(2);
  1818. !%%%%%%%%%%%%%%%%%%%%%%
  1819. ! Cell definition
  1820. !%%%%%%%%%%%%%%%%%%%%%%
  1821. ! eta(1) : prey
  1822. ! eta(2) : predator
  1823. f_set deta_tef(1) = apar*eta(1)-apar*ff(1);
  1824. f_set deta_tef(2) = - cpar*eta(2) + cpar*ff(1);
  1825. @end example
  1826. The starting points for cells are entered like:
  1827. @example
  1828. ! initial state
  1829. ! -------------
  1830. eta(1) = 1.;
  1831. eta(2) = 1.;
  1832. @end example
  1833. If there are observations, they are entered as special transferts with
  1834. index above @code{mp}, for example:
  1835. @example
  1836. f_set Phi_tef(mp+1) = ff(1) ;
  1837. @end example
  1838. @node Programming with cmz directives
  1839. @section Programming with cmz directives
  1840. @menu
  1841. * Cmz directives used with @Minik{}::
  1842. * Using cmz directives in @Minik{}::
  1843. @end menu
  1844. @node Cmz directives used with @Minik{}
  1845. @subsection Cmz directives used with @Minik{}
  1846. The main feature of cmz directive is to use code conditionnaly for a given
  1847. select flag. For example when the double precision is selected
  1848. (@pxref{Double precision}) the use of the conditionnal
  1849. @code{double} flag may be required in case there is a different subroutine
  1850. name for different types. If, for example, the user use the subroutine
  1851. @code{smysub} for simple precision and @code{dmysub} for double
  1852. precision the following code is an example of what could appear in the
  1853. user code:
  1854. @verbatim
  1855. +IF,double
  1856. call dmysub(eta);
  1857. +ELSE
  1858. call smysub(eta);
  1859. +ENDIF
  1860. @end verbatim
  1861. For a complete reference on cmz directives see the appendix
  1862. @ref{Cmz directives reference}.
  1863. @node Using cmz directives in @Minik{}
  1864. @subsection Using cmz directives in @Minik{}
  1865. In cmz the KEEP and DECK have their cmz directives preprocessed as part
  1866. of the source files extraction. And the +KEEP and +DECK
  1867. directives are automatically
  1868. set when creating the KEEP or DECK. With make, files with these directives
  1869. has to be created within the files that are to be preprocessed by the
  1870. cmz directives preprocessor.
  1871. To be processed by make, a file that contains cmz directives
  1872. should have a file suffix corresponding
  1873. with the language of the resulting file and with the normal file suffix of
  1874. that language. More precisely @samp{cm} should be added before the normal
  1875. file suffix and after the @samp{.}. Therefore if the resulting file language
  1876. is associated with a suffix @samp{.@var{suf}}, the file with cmz directives
  1877. should have a @samp{.cm@var{suf}} suffix. The tradition is to have
  1878. a different suffix for main files and include files.
  1879. To add directories searched for @dfn{cmfiles} (files with cmz directives)
  1880. they should be added to the @code{CMFDIRS} makefile variable, separated
  1881. by @samp{:}.
  1882. Rules for preprocessing of the files are defined in the file
  1883. @file{Makefile.miniker} for the file types described in
  1884. @ref{tab:cmfile_suffix}:
  1885. @float table, tab:cmfile_suffix
  1886. @multitable {fortran preprocessed} {include/keep} {cmfile suffix} {suffix} {language}
  1887. @headitem language @tab file type @tab cmfile suffix @tab suffix @tab language
  1888. @item fortran @tab main/deck @tab .cmf @tab .f @tab ftn
  1889. @item fortran preprocessed @tab main/deck @tab .cmF @tab .F @tab f77
  1890. @item fortran preprocessed @tab include/keep @tab .cminc @tab .inc @tab f77
  1891. @item mortran @tab main/deck @tab .cmmtn @tab .mtn @tab mtn
  1892. @item mortran @tab include/keep @tab .cmmti @tab .mti @tab mtn
  1893. @end multitable
  1894. @caption{Association between file language, file type, file suffixes and
  1895. language identifier in cmz directives. A main file is called a @dfn{deck}
  1896. in cmz and an include file is called a @dfn{keep}.}
  1897. @end float
  1898. @node Dynamic system analysis
  1899. @chapter Dynamic analysis of systems in @Minik{}
  1900. @menu
  1901. * Sensitivities::
  1902. * Adjoint model and optimisation::
  1903. * Kalman filter::
  1904. * Feedback gain::
  1905. * Stability of fastest modes::
  1906. * Generalized TLS::
  1907. @end menu
  1908. @node Sensitivities
  1909. @section Automatic sensitivity computation
  1910. @cindex sensitivities
  1911. An obvious advantage of having acces to the Jacobian matrices along the
  1912. system trajectory concerns automatic sensitivity analyses, as either:
  1913. @itemize @bullet
  1914. @item the sensitivity of all variables to perturbation in the initial condition
  1915. of one state variable;
  1916. @item the same sensitivities to an initial pulse (or step) on a transfer;
  1917. @item the same sensitivities to a series of pulses (or steps) on a transfer;
  1918. @item the same for a change in a parameter, eventually during the run;
  1919. @item the sensitivity of the matrix of advance in state space to a change
  1920. in a parameter.
  1921. @end itemize
  1922. This is declared in Zinit as:
  1923. @example
  1924. ! -------------
  1925. ! Sensitivities
  1926. ! -------------
  1927. Sensy_to_var
  1928. < var: eta_pray, pert: INIT;
  1929. var: eta_pred, pert: INIT;
  1930. >;
  1931. @end example
  1932. Each variable at origin of a perturbation is declared as @code{var:},
  1933. and the type of perturbation in @code{pert:}. Here, INIT conditions are
  1934. only allowed because the two variables are states variables. For transfers,
  1935. @code{pert: pulse} corresponds to an initial pulse, @code{pert: step_resp}
  1936. and @code{pert: step_eff} to initial steps, the difference between
  1937. @code{_resp} (response form)
  1938. and @code{_eff} (effect form) concerns the
  1939. diagonal only of the sensitivity matrix
  1940. (see Feedback gains in non-linear models).
  1941. Non initial perturbation can also be asked for:
  1942. @example
  1943. Sensy_to_var
  1944. <
  1945. !* var: eta_courant_L, pert: init at 100;
  1946. !* var: ff_T_czcx, pert: pulse at 100 every 20;
  1947. !* var: ff_Psi_Tczcx, pert: step_eff;
  1948. !* var: ff_Psi_Tczcx, pert: step_Resp at 10 every 100;
  1949. ! *** premiers tests identiques a lorhcl.ref
  1950. var: ff_courant_L , pert: step_eff;
  1951. var: ff_T_czcx , pert: step_eff;
  1952. var: ff_Psi_Tczcx , pert: step_eff;
  1953. var: ff_Psi_Tsz , pert: pulse at 100 every 50;
  1954. >;
  1955. @end example
  1956. In this example taken from @file{lorhcl}, a sensitivity can increase so as to
  1957. trespass the Fortran capacity, so that each sensitivity vector (matrix column)
  1958. can be reset at some time-increment @code{at III every JJJ;}
  1959. It is noteworthy that these sensitivity analyses are not based
  1960. on difference between two runs with different initial states or
  1961. parameter values, but on the formal derivatives of the model. This method
  1962. is not only numerically robust, but is also rigorously funded as based on
  1963. the TLS of the model@footnote{For a short introduction to automatic
  1964. sensitivity analysis, see the document:@*
  1965. @url{http://lmd.jussieu.fr/zoom/doc/sensibilite.ps}, in French,
  1966. or ask for the more complete research document to a member of the TEF-ZOOM
  1967. collaboration}.
  1968. If the @code{dimetaphi} sequence is built by the users, he should declare
  1969. the number of perturbing variables as @code{nxp=}:
  1970. @example
  1971. parameter (nxp=np,nyp=0,nzp=0);
  1972. @end example
  1973. here, all state variables are considered as perturbing variables.
  1974. @cindex sensitivity, output
  1975. @cindex output, sensitivity
  1976. @cindex @file{sens.data}
  1977. @cindex @file{sigma.data}
  1978. The sensitivity vectors are output in the result files @file{sens.data} for
  1979. cells and @file{sigma.data} for transfers. In those files the first column
  1980. corresponds again with time, and the other columns are relative sensitivities of the cell
  1981. states (in @file{sens.data}) and transfers (in @file{sigma.data})
  1982. with respect to the initial value of the perturbed state.
  1983. In our predator-prey example, the second column of @file{sens.data} will contain
  1984. the derivative of @math{\eta_1(t)} with respect to @math{\eta_1(t=0)}.
  1985. Drawing the
  1986. second column of @file{sens.data} against the first one
  1987. gives the time evolution of the sensitivity of @code{eta-pred}
  1988. to a change in the initial value of @code{eta-pray}. One can check
  1989. in that it is set to 1 at @math{t=0}:
  1990. @example
  1991. # Sensy_to: eta_pray 3 eta_pred 5
  1992. # time \\ of: eta_pray eta_pred eta_pray eta_pred
  1993. 0.00000E+00 1.00000E+00 0.00000E+00 0.00000E+00 1.00000E+00
  1994. 1.00000E-02 9.90868E-01 1.11905E-02 -1.26414E-02 9.98859E-01
  1995. @end example
  1996. The two last columns are the state sensitivity to a change in initial conditions
  1997. of the number of predators.
  1998. In the same way, the @var{j+1}th column of @file{sigma.data} will be the
  1999. derivative of @math{\phi_{j}(t)} with respect to @math{\eta_i(t=0)}. Here:
  2000. @example
  2001. # Sensy_to: eta_pray eta_pred
  2002. # time \\ of: ff_interact ff_interact
  2003. 0.00000E+00 1.60683E+00 8.47076E-01
  2004. 1.00000E-02 1.59980E+00 8.18164E-01
  2005. @end example
  2006. the unique transfer variable gives rise to two sensitivity columns.
  2007. Sensitivity studies are usefull to assess the
  2008. predictability properties of the corresponding system.
  2009. @menu
  2010. @c * Initial state sensitivity::
  2011. @c * Sensitivity to a pulse or a step on transfer::
  2012. @c * Extended Sensitivity studies::
  2013. * Sensitivity to a parameter::
  2014. * Advance matrix sensitivity::
  2015. @end menu
  2016. @node Sensitivity to a parameter
  2017. @subsection Sensitivity to a parameter
  2018. A forward sensitivity to a parameter will be computed when specified as
  2019. described in @ref{Parameters}. For example, suppose that
  2020. the sensitivity to an initial change in the @code{apar} parameter of
  2021. the predator model is of interest.
  2022. @c In that case the number of
  2023. @c parameters should be set to 1 in @file{dimetaphi}:
  2024. @c
  2025. @c @example
  2026. @c parameter (lp=2);
  2027. @c @end example
  2028. The sensitivity calculs is turned on as a forward
  2029. parameter specified on the @code{Free_parameter} list:
  2030. @example
  2031. Free_parameter: [fwd: apar, cpar];
  2032. @end example
  2033. The result are in @file{sensp.data} for cells and @file{sigmap.data}
  2034. for transfers.
  2035. @example
  2036. # Sensy_to: pi_prandtl 3 4 pi_rayleigh_ 6
  2037. # time \\ of: eta_courant_ eta_T_czcx eta_T_sz eta_courant_ eta_T
  2038. 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.000
  2039. 2.00000E-03 -4.77172E-03 -3.99170E-05 3.55971E-05 -9.94770E-05 -1.004
  2040. @end example
  2041. In the above example from @file{lorhcl} sensitivity of the three states with respect
  2042. to an initial change in two parameters are independantly given (first line also numbers
  2043. the column to easy gnuplot using).
  2044. @node Advance matrix sensitivity
  2045. @subsection Advance matrix sensitivity
  2046. It is possible to look at the sensitivity of the matrix of advance in
  2047. states space (the matrix @code{aspha}) with regard to a parameter.
  2048. The parameter must be accounted for in the parameter number and be in the
  2049. parameter list, flagged as the matrix @code{mx} parameter, like in
  2050. @example
  2051. Free_parameter: [mx: apar], cpar;
  2052. @end example
  2053. @vindex d_pi_aspha(.,.)
  2054. This feature is associated with a selecting flag, @samp{dPi_aspha}. One gets
  2055. the result in the matrix @code{d_pi_aspha(.,.)} of dimension
  2056. (@code{np},@code{np}).
  2057. This matrix may be used to compute other quantities, for example
  2058. it may be used to compute the sensitivity of the eigenvalues of
  2059. the state-advance matrix with regard to the @code{[fwd]} parameter.
  2060. These additional computations have to be programmed by the user in
  2061. @file{zsteer} with matrices declared and initialized in
  2062. @file{zinit}. An example is given in the example @file{lorhcl}
  2063. provided with the @Minik{} installation files, following a method proposed
  2064. by Stephane Blanco.
  2065. @node Adjoint model and optimisation
  2066. @section Adjoint model and optimisation with @Minik{}
  2067. In the following a possible use of @Minik{} for optimisation is discussed.
  2068. More precisely the use of adjoint and control laws in @Minik{} are presented.
  2069. Optimisation isn't the only application of these tools, but it is the most
  2070. common one. In that case the adjoint may be used to determine the gradient of a
  2071. functional to perturbations in the control laws, and an optimisation process
  2072. can use this
  2073. information to search for the optimum.
  2074. Another application of the adjoint is to compute the sensitivity of a
  2075. cost function to parameters (the ones declared in the @code{free_parameters:}' list.
  2076. Note that the cost function can be sensitive to probe's variables, even if these are
  2077. uncoupled with standard variables in the forward calculations; this is the case
  2078. when minimizing a quadratic distance function between probes (from the model)
  2079. and the corresponding measurements.
  2080. The code is close transcription of the mathematical calculus described
  2081. in@* @url{http://www.lmd.jussieu.fr/ZOOM/doc/Adjoint.pdf} . It essentialy reverse time and
  2082. transpose the four Jacobian matrices: states and transfers are saved in array dimensionned
  2083. with @code{maxstep} Fortran parameter.
  2084. @menu
  2085. * Overview of optimisation with @Minik{}::
  2086. * Control laws::
  2087. * Cost function coding and adjoint modeling::
  2088. * Sensitivity of cost function to parameters::
  2089. @end menu
  2090. @node Overview of optimisation with @Minik{}
  2091. @subsection Overview of optimisation with @Minik{}
  2092. @cindex adjoint
  2093. @cindex optimisation
  2094. In the proposed method, @Minik{} is run twice, one time forward and then
  2095. backward to determine the trajectory and the adjoint model. After that the
  2096. control laws are modified by a program external to @Minik{}. The same steps
  2097. are repeated until convergence. More pecisely,
  2098. @table @strong
  2099. @item forward
  2100. The command law @math{h(t)} is given (by an explicit law or taken from a file).
  2101. The trajectory is computed in a classical way, with the additionnal computation
  2102. of the functional to be optimised, @math{J}, prescribed with specific
  2103. @code{f_set} macros. The states, transfers and control laws are stored.
  2104. @item backward
  2105. The adjoint variable is computed from the last time @math{T} backward. The
  2106. time increment is re-read as it could have changed during the forward
  2107. simulation. The system is solved by using the same technics as in the forward
  2108. simulation, but with a negative time step.
  2109. @item external phase
  2110. Now the command should be corrected. This step isn't covered here, but, for
  2111. example, minuit the optimisation tool from the CERN could be used.
  2112. In order to ease such a use of @Minik{}, the principal program has to be
  2113. compiled as a subroutine to be driven by an external program
  2114. (@pxref{Calling the model code}).
  2115. @end table
  2116. The functionnal @math{J} to be optimised is defined as
  2117. @ifset latex
  2118. @tex
  2119. \begin{equation}
  2120. J = \psi[\eta(T),\varphi(T) ,h(T)] + \int_0 ^T {l[\eta(\tau),\varphi(\tau),h(\tau)]}\, d\tau
  2121. \end{equation}
  2122. @end tex
  2123. @end ifset
  2124. @ifclear latex
  2125. @tex
  2126. $$
  2127. J = \psi[\eta(T),\varphi(T) ,h(T)] + \int_0 ^T {l[\eta(\tau),\varphi(\tau),h(\tau)]}\, d\tau
  2128. $$
  2129. @end tex
  2130. @ifnottex
  2131. @noindent @math{J = psi(eta(T),phi(T),h(T)) + int_0^T l(eta(tau),phi(tau),h(tau)) d tau}
  2132. @end ifnottex
  2133. @end ifclear
  2134. @cindex final cost
  2135. @cindex integrand cost
  2136. Where @math{\psi} is the final cost function, @math{l} is the integrand
  2137. cost function and @math{h} represents the control laws variations.
  2138. The general use of the adjoint model of a system is the determination of the
  2139. gradient of this @math{J} functional to be optimised, with respect to perturbations
  2140. of the original conditions of the reference trajectory, that is, along its
  2141. GTLS@footnote{General Tangent Linear System, i.e. the TLS circulating along a trajectory.
  2142. See the explanation in the document
  2143. @url{http://www.lmd.jussieu.fr/Zoom/doc/Adjoint.pdf} (in French).}.
  2144. @node Control laws
  2145. @subsection Control laws
  2146. @vindex zcommand
  2147. @cindex command law
  2148. Each control law is associated with one cell or transfer equation, meaning that a command
  2149. associated with an equation does not appear in any other equation.
  2150. It is still possible
  2151. to add commands acting anywhere by defining a transfer equal to that command.
  2152. The control laws associated with states are in the @code{ux_com(.)} array,
  2153. control laws associated with transfers are in the @code{uy_com(.)} array.
  2154. The control laws may be prescribed even when there is no adjoint computed,
  2155. nor any optimisation, and they are used during simulation, in which case they will
  2156. act as external sources. To enable
  2157. the use of commands, the logical flag @code{Zcommand} should be @code{.true.}.
  2158. @cindex @file{uxcom.data}
  2159. @cindex @file{uycom.data}
  2160. The command can be given either as:
  2161. @enumerate
  2162. @item a table of numerical
  2163. values in the files @file{uxcom.data} and @file{uycom.data}.
  2164. @item a function
  2165. @vindex zlaw
  2166. @cindex @file{zcmd_law}
  2167. @cindex @file{zcmd_law.inc}
  2168. of the problem variables. To turn that feature on the logical flag
  2169. @code{Zlaw} should be set to @code{.true.} in @file{zinit}. The sequence
  2170. @file{zcmd_law} should hold
  2171. the code filling the @code{ux_com(.)} and @code{uy_com(.)} arrays, as the code
  2172. from that sequence is used whenever the control laws are needed.
  2173. In that case the files @file{uxcom.data} and @file{uycom.data} will
  2174. be filled by the command values generated by the function along the trajectory.
  2175. @end enumerate
  2176. For example in the Lotka-Volterra model, the parameter @code{apar} could
  2177. be a control variable.
  2178. In that case, @code{apar} would be defined as the variable @code{ux_com(1)},
  2179. and either entered as a law
  2180. in the sequence @file{zcmd_law} , either written in the file @file{uxcom.data}
  2181. step by step. In that case, there must be a perfect corresponodence between time
  2182. of the commands and time of the run.
  2183. @node Cost function coding and adjoint modeling
  2184. @subsection Cost function coding and adjoint modeling
  2185. @vindex zback
  2186. @findex cout_Psi
  2187. @findex cout_l
  2188. First of all the flag @code{zback} should be set to @code{.true.} in order to
  2189. allow adjoint model computation:
  2190. @example
  2191. Zback=.true.;
  2192. @end example
  2193. The two functions @code{cout_Psi} corresponding with the final cost and
  2194. @code{cout_l} corresponding with the integrand cost are set up with the
  2195. @code{f_set} macros.
  2196. @defmac f_set cout_Psi = f(eta(.),ff(.),ux_com(.),uy_com(.))
  2197. This macro defines the final cost function.
  2198. @code{f} is a fortran
  2199. expression which may be function of cell state variables,
  2200. @samp{eta(1)}@dots{}@samp{eta(np)}, transfers
  2201. @samp{ff(1)}@dots{}@samp{ff(mp)},
  2202. state control laws
  2203. @samp{ux_com(1)}@dots{}@samp{ux_com(np)}, and transfer control laws
  2204. @samp{uy_com(1)}@dots{}@samp{uy_com(mp)}.
  2205. @end defmac
  2206. @defmac f_set cout_l = f(eta(.),ff(.),ux_com(.),uy_com(.))
  2207. This macro defines the integrand cost function.
  2208. @code{f} is a fortran
  2209. expression which may be function of cell state variables,
  2210. @samp{eta(1)}@dots{}@samp{eta(np)}, transfers
  2211. @samp{ff(1)}@dots{}@samp{ff(mp)},
  2212. state control laws
  2213. @samp{ux_com(1)}@dots{}@samp{ux_com(np)}, and transfer control laws
  2214. @samp{uy_com(1)}@dots{}@samp{uy_com(mp)}.
  2215. @end defmac
  2216. For example, the following code sets a cost function for the masselottes
  2217. model:
  2218. @example
  2219. ! Initialisation
  2220. F_set cout_Psi = eta_move(inode:1);
  2221. !and f_set cout_l integrand in the functionnal
  2222. F_set cout_l = 0.;
  2223. @end example
  2224. In that example the functional is reduced to the final value
  2225. of the first state component.
  2226. Here, the adjoint vector will correspond to the final sensitivity
  2227. (at @math{t=0}) of
  2228. that component (here the first masselotte position) to a perturbation in
  2229. all initial conditions@footnote{For detailed explanation of the adjoint model,
  2230. see the document in
  2231. @uref{http://www.lmd.jussieu.fr/@/ZOOM/doc/Adjoint.pdf,pdf}
  2232. or @uref{http://www.lmd.jussieu.fr/@/ZOOM/doc/Adjoint.pdf,.ps.gz}}.
  2233. @c In the code, the variables @code{v_adj(.)} and @code{w_adj(.)}
  2234. @c are respectively adjoint to @code{eta(.)} and @code{ff(.)}. They are written in the
  2235. @c two files: @file{vadj.data} and @file{wadj.data}.
  2236. The following variables are set during the backward phase, and output
  2237. in the associated files:
  2238. @multitable {@code{gradufj(.)}} {@file{hamilton.data}} {time increment, hamiltonian, cost function increment}
  2239. @headitem var @tab file @tab explanation
  2240. @c @item @code{} @tab @file{.data} @tab @tab
  2241. @item @code{v_adj(.)} @tab @file{vadj.data} @tab adjoint to @code{eta(.)}
  2242. @item @code{w_adj(.)} @tab @file{wadj.data} @tab adjoint to @code{ff(.)}
  2243. @item @code{wadj(mp+.)} @tab @file{gradmuj.data} @tab adjoint to @code{ff(mp+.)}
  2244. @item @code{graduej(.)} @tab @file{gradxj.data} @tab adjoint to @code{ux_com(.)}
  2245. @item @code{gradufj(.)} @tab @file{gradyj.data} @tab adjoint to @code{uy_com(.)}
  2246. @item @code{hamilton} @tab @file{hamilton.data} @tab time increment, hamiltonian, cost function increment
  2247. @end multitable
  2248. @node Sensitivity of cost function to parameters
  2249. @subsection Sensitivity of cost function to parameters
  2250. @cindex @file{gradpj.data}
  2251. The sensitivity of the cost function to all the parameters given as
  2252. arguments of @code{Free_parameters} is computed. For the
  2253. predator model the sensitivity of a cost function consisting in
  2254. the integral of the predator population with respect with
  2255. @code{apar} an @code{cpar} is obtained with a number of parameters
  2256. set to 2 in @file{dimetaphi}:
  2257. @example
  2258. parameter (lp=2);
  2259. @end example
  2260. And the cost function and @code{Free_parameters} list in @file{zinit}:
  2261. @example
  2262. f_set cout_Psi = eta(2);
  2263. f_set cout_l = eta(2);
  2264. Free_parameters: apar,cpar;
  2265. @end example
  2266. @code{apar} and @code{cpar} also have to be given a value.
  2267. The result is output in @file{gradpj.data}.
  2268. @node Kalman filter
  2269. @section Kalman filter
  2270. @cindex Kalman filter
  2271. @cindex variance-covariance matrices, general
  2272. @cindex observations, general
  2273. The Kalman filter allows for data assimilation along the model run. In
  2274. that case it is assumed that there is a real-world model with stochastic
  2275. perturbations on the states, and that noisy observations are available.
  2276. The situation implemented in @Minik{} corresponds to a continuous
  2277. stochastic perturbation on the state, and discrete noisy observations.
  2278. In the @acronym{TEF} this leads to:
  2279. @ifset latex
  2280. @tex
  2281. \begin{eqnarray}
  2282. \partial_t \eta (t) &=& g(\eta(t),\varphi(t)) + W(t) \mu\\
  2283. \varphi(t) &=& f(\eta(t),\varphi(t))\\
  2284. \omega(t) &=& h ( \eta(t) , \varphi(t)) + \nu
  2285. \end{eqnarray}
  2286. @end tex
  2287. @end ifset
  2288. @ifclear latex
  2289. @tex
  2290. $$\eqalign{
  2291. \partial_t \eta (t) &= g(\eta(t),\varphi(t)) + W(t) \mu\cr
  2292. \varphi(t) &= f(\eta(t),\varphi(t))\cr
  2293. \omega(t) &= h ( \eta(t) , \varphi(t)) + \nu\cr
  2294. }$$
  2295. @end tex
  2296. @ifnottex
  2297. @noindent @math{d eta(t)/d t = g(eta(t),phi(t)) + W(t) mu@*
  2298. phi(i) = f(eta(t),phi(t))@*
  2299. omega(t) = h(eta(t), phi(t)) + nu }
  2300. @end ifnottex
  2301. @end ifclear
  2302. @c FIXME partout omega
  2303. @c (notice that in this paragraph, $\omega$ stands for the probe vector $\mu$ elsewhere,
  2304. @c and $\mu$ is here a noise source.
  2305. The observations @math{\omega} are available at discrete time steps @math{t=s_i}. The
  2306. stochastic perturbation on state, @math{\mu} is characterized by a
  2307. variance-covariance matrix @math{Q} and the noise on the observation,
  2308. @math{\nu} has a variance-covariance matrix @math{R}. @math{W} relates states
  2309. with stochastic perturbations. At each time step the Kalman filter recomputes
  2310. an estimation of the state and the variance-covariance matrix of the state.
  2311. In the following we use the example of a linear model with perturbation
  2312. on state and observation of state. The model has 3 states and 3 corresponding
  2313. transfers (equal to the states), but the error on the state is of dimension
  2314. 2. The 3 states are observed. The corresponding equations read:
  2315. @ifset latex
  2316. @tex
  2317. \begin{equation}
  2318. \left\{
  2319. \begin{array}{cc}
  2320. \partial_t \eta_1 =& a_{11} \eta_1 + a_{12} \varphi_2 + a_{13} \varphi_3 + W_{11} \mu_1 + W_{12} \mu_2\\
  2321. \partial_t \eta_2 =& a_{21} \varphi_1 + a_{22} \eta_2 + a_{23} \varphi_3 + W_{21} \mu_1 + W_{22} \mu_2\\
  2322. \partial_t \eta_3 =& a_{31} \varphi_1 + a_{32} \varphi_2 + a_{33} \eta_3 + W_{31} \mu_1 + W_{32} \mu_2
  2323. \end{array}
  2324. \right.
  2325. \end{equation}
  2326. \begin{equation}
  2327. \left\{
  2328. \begin{array}{cc}
  2329. \varphi _1 =& \eta _1\\
  2330. \varphi _2 =& \eta _2\\
  2331. \varphi _3 =& \eta _3
  2332. \end{array}
  2333. \right.
  2334. \end{equation}
  2335. \begin{equation}
  2336. \left\{
  2337. \begin{array}{cc}
  2338. \omega _1 =& \varphi _1 + \nu_1\\
  2339. \omega _2 =& \eta _2 + \nu_2 \\
  2340. \omega _3 =& \eta _3 + \nu_3
  2341. \end{array}
  2342. \right.
  2343. \end{equation}
  2344. @end tex
  2345. @end ifset
  2346. @ifclear latex
  2347. @tex
  2348. $$\left\{\eqalign{
  2349. \partial_t \eta_1 &= a_{11} \eta_1 + a_{12} \varphi_2 + a_{13} \varphi_3 + W_{11} \mu_1 + W_{12} \mu_2\cr
  2350. \partial_t \eta_2 &= a_{21} \varphi_1 + a_{22} \eta_2 + a_{23} \varphi_3 + W_{21} \mu_1 + W_{22} \mu_2\cr
  2351. \partial_t \eta_3 &= a_{31} \varphi_1 + a_{32} \varphi_2 + a_{33} \eta_3 + W_{31} \mu_1 + W_{32} \mu_2
  2352. }\right.$$
  2353. $$\left\{\eqalign{
  2354. \varphi _1 &= \eta _1\cr
  2355. \varphi _2 &= \eta _2\cr
  2356. \varphi _3 &= \eta _3
  2357. }\right.$$
  2358. $$\left\{\eqalign{
  2359. \omega _1 &= \varphi _1 + \nu_1\cr
  2360. \omega _2 &= \eta _2 + \nu_2 \cr
  2361. \omega _3 &= \eta _3 + \nu_3
  2362. }\right.$$
  2363. @end tex
  2364. @ifnottex
  2365. Cells:@*
  2366. @noindent @math{d eta_1/dt = a_11 eta_1 + a_12 phi_2 + a_13 phi_3 + W_11 mu_1 + W_12 mu_2@*
  2367. d eta_2/dt = a_21 phi_1 + a_22 eta_2 + a_23 phi_3 + W_21 mu_1 + W_22 mu_2@*
  2368. d eta_3/dt = a_31 phi_1 + a_32 phi_2 + a_33 eta_3 + W_31 mu_1 + W_32 mu_2}
  2369. Transfers:@*
  2370. @noindent @math{phi_1 = eta_1@*
  2371. phi_2 = eta_2@*
  2372. phi_3 = eta_3}
  2373. Observations:@*
  2374. @noindent @math{omega_1 = phi_1 + nu_1@*
  2375. omega_2 = eta_2 + nu_2@*
  2376. omega_3 = eta_3 + nu_3}
  2377. @end ifnottex
  2378. @end ifclear
  2379. @menu
  2380. * Coding the Kalman filter::
  2381. * Kalman filter run and output::
  2382. * Executing code after the analysis::
  2383. @end menu
  2384. @node Coding the Kalman filter
  2385. @subsection Coding the Kalman filter
  2386. @vindex zkalman
  2387. First of all the Kalman filter code should be activated. The observations
  2388. code is also required (@pxref{Observations}).
  2389. If cmz is used the code
  2390. should be selected with the select flag kalman
  2391. in the @file{selseq.kumac}:
  2392. @example
  2393. sel kalman
  2394. @end example
  2395. With make the @code{kalman} variable should be set to 1:
  2396. @example
  2397. kalman = 1
  2398. @end example
  2399. The kalman code is actually used by setting the flag
  2400. @code{zkalman} to @code{.true.}, for example in the @file{zinit}:
  2401. @example
  2402. zkalman = .True.;
  2403. @end example
  2404. @c This will set the @code{zobs} and @code{zdata} flags to @code{.true.}
  2405. @c (@pxref{Observations and data}).
  2406. With the Kalman filter the dimension of estimated states, of the error
  2407. on the state and of the
  2408. observation, the @math{W} matrix, the observation function,
  2409. the initial
  2410. variance-covariance matrices on the state and the variance-covariance matrices
  2411. of errors have to be given.
  2412. @menu
  2413. * Kalman filter vectors dimensions::
  2414. * Error and observation matrices::
  2415. @end menu
  2416. @node Kalman filter vectors dimensions
  2417. @subsubsection Kalman filter vectors dimensions
  2418. @cindex error vector dimension
  2419. @cindex @file{dimetaphi}, Kalman filter
  2420. These dimensions should be set in the @file{zinit} sequence.
  2421. The size of the estimated states is given by the parameter @code{nkp}.
  2422. You can set this to @code{np} if all the states are estimated, but in case
  2423. there are some deterministic state variables, @code{nkp} may be less than
  2424. @code{np}. In that case the first @code{nkp} elements of @code{eta(.)}
  2425. will be estimated using the Kalman filter.
  2426. The error on state dimension is associated with the parameter @code{nerrp}
  2427. and the size of the observations vector is @code{mobs}
  2428. (@pxref{Observations}). In our example the dimensions are set with:
  2429. @example
  2430. parameter (nkp=np);
  2431. parameter (mobs=3);
  2432. parameter (nerrp=2);
  2433. @end example
  2434. All the states are estimated,
  2435. there are 3 observation functions and the error on the state vector is of
  2436. dimension 2.
  2437. If the sizes are set explicitely, the parameters should be set in
  2438. @file{dimetaphi}.
  2439. @node Error and observation matrices
  2440. @subsubsection Error and observation matrices
  2441. @cindex variance-covariance matrices
  2442. @cindex observations
  2443. @cindex @file{zinit}, Kalman filter
  2444. @subsubheading Initial variance-covariance matrix on the state
  2445. @cindex initial variance-covariance on states
  2446. @vindex covfor(.,.)
  2447. The variance-covariance on the state matrix is @code{covfor(.,.)}. The initial
  2448. values have to be given for this matrix, as in our example:
  2449. @example
  2450. covfor(1,1) = 1000.; covfor(1,2) = 10.; covfor(1,3) = 10.;
  2451. covfor(2,1) = 10.; covfor(2,2) = 5000.; covfor(2,3) = 5.;
  2452. covfor(3,1) = 10.; covfor(3,2) = 5.; covfor(3,3) = 2000.;
  2453. @end example
  2454. This matrix is updated by the filter at each time step because the states
  2455. are pertubated by some noise, and when assimilation takes place as new
  2456. information reduce the error.
  2457. @subsubheading Observations and error on state matrix
  2458. @cindex variance-covariance matrix on state
  2459. @vindex mereta(.,.)
  2460. The matrix that relates errors on states vector components to states,
  2461. corresponding with @math{W} is @code{mereta(.,.)}. In our example it is
  2462. set by:
  2463. @example
  2464. mereta(1,1) = 1.; mereta(1,2) = 0.;
  2465. mereta(2,1) = 0.; mereta(2,2) = 1.;
  2466. mereta(3,1) = 0.5; mereta(3,2) = 0.5;
  2467. @end example
  2468. The observation functions are set by a @code{f_set} macro with
  2469. @code{Obs_tef(.)} as described in @ref{Observations}.
  2470. In our example the observation functions are set by:
  2471. @example
  2472. f_set Obs_tef(1) = ff(1) ;
  2473. f_set Obs_tef(2) = eta(2);
  2474. f_set Obs_tef(3) = eta(3);
  2475. @end example
  2476. @subsubheading Error variance-covariance matrices
  2477. @cindex variance-covariance error
  2478. @vindex covobs(.,.)
  2479. The variance-covariance matrix on observation noise is @code{covobs(.,.)}
  2480. set, in our example, by:
  2481. @example
  2482. covobs(1,1) = 0.3; covobs(1,2) = 0.; covobs(1,3) = 0.;
  2483. covobs(2,1) = 0.; covobs(2,2) = 0.1; covobs(2,3) = 0.;
  2484. covobs(3,1) = 0.; covobs(3,2) = 0.; covobs(3,3) = 0.2;
  2485. @end example
  2486. @vindex coveta(.,.)
  2487. The variance-covariance matrix on state noise is @code{coveta(.,.)}
  2488. set, in our example, by:
  2489. @example
  2490. coveta(1,1) = 0.2; coveta(1,2) = 0.001;
  2491. coveta(2,1) = 0.001; coveta(2,2) = 0.1;
  2492. @end example
  2493. These matrices are not changed during the run of the model as part
  2494. of the filtering process. They may be changed by the user in @file{zsteer}.
  2495. @node Kalman filter run and output
  2496. @subsection Kalman filter run and output
  2497. @menu
  2498. * Feeding the observations::
  2499. * Kalman filter results::
  2500. @end menu
  2501. @node Feeding the observations
  2502. @subsubsection Feeding the observations to the model
  2503. @vindex vobs(.)
  2504. @vindex zgetobs
  2505. @cindex @file{zsteer}, Kalman filter
  2506. The observations must be made available to the model during the run. These
  2507. observations are set in the @code{vobs(.)} array, and the assimilation
  2508. (also called the analysis step of the filter) takes
  2509. place if the logical variable @code{zgetobs} is @code{.true.}
  2510. (@pxref{Data}).
  2511. These steps are
  2512. typically performed in the @file{zsteer} sequence. In this sequence there should
  2513. be some code such that when there are data ready to
  2514. be assimilated, @code{zgetobs} is set to @code{.true.} and the data is
  2515. stored in @code{vobs(.)}, ready for the next step processing.
  2516. @node Kalman filter results
  2517. @subsubsection Kalman filter results
  2518. @cindex results, Kalman filter
  2519. @cindex Kalman filter results
  2520. @cindex output, Kalman filter
  2521. @cindex Kalman filter output
  2522. @cindex @file{data.data}
  2523. The estimated states and transfers are still in the same @samp{.data} files,
  2524. @file{res.data} and @file{tr.data} and there is the additional file with
  2525. observations, called @file{obs.data} (@pxref{Observations}).
  2526. Each time @code{zgetobs} is @code{.true.} the data, and the optimally
  2527. weighted innovations are output
  2528. in the file associated with data, @file{data.data} (@pxref{Data}).
  2529. @node Executing code after the analysis
  2530. @subsection Executing code after the analysis
  2531. The analysis takes place before the time step advance when @code{zgetobs}
  2532. is @code{.true.}. It may be usefull to add some code after the analysis
  2533. and before the time step advance. For example the analysis may lead to
  2534. absurd values for some states or parameters, it could be usefull to correct
  2535. them in that case. The sequence included after the analysis is called
  2536. @file{kalsteer}. At this point, in addition to the usual variables
  2537. the following variables could be usefull:
  2538. @vtable @code
  2539. @item etafor(.)
  2540. The state before the analysis.
  2541. @item kgain(.)
  2542. The Kalman gain.
  2543. @item innobs(.)
  2544. The innovation vector (observations coherent with the states minus data
  2545. values).
  2546. @item covana(.,.)
  2547. The variance-covariance error matrix after the analysis.
  2548. @end vtable
  2549. At each time step the derivative of the observation function with respect
  2550. to transfer and cells variables are recomputed. The elimination of
  2551. transfers is also performed to get the partial derivative of the observation
  2552. function of the equivalent model, with states only, with respect to the
  2553. states. In other words, the Kalman filter does not follow the TEF formalism, because
  2554. the advance of the var-covar matrix could not yet be set in the TEF form.
  2555. @c There is a corresponding additional matrix:
  2556. @vtable @code
  2557. @c @item obetad(.,.)
  2558. @c derivative of observation function with respect to transfers.
  2559. @c @item obphid(.,.)
  2560. @c derivative of observation function with respect to cell variables.
  2561. @item obspha(.,.)
  2562. derivative of observation function in state space with respect to
  2563. cell variables.
  2564. @end vtable
  2565. @node Feedback gain
  2566. @section Feedback gain
  2567. @cindex Borel sweep
  2568. @cindex Feedback gain
  2569. The feedback dynamic gain associated with a feedback loop
  2570. can be expressed as the inverse Borel
  2571. transform of the coefficient of the reduced scalar
  2572. coupling matrix, @math{g(\tau)},
  2573. associated with a transfer.
  2574. A Borel sweep provides this @math{g(\tau)}. Therefore it is
  2575. an interesting tool for the characterization of the feedback loop@footnote{
  2576. More generally, the Borel sweep allows
  2577. the numerical study of the dependency in @math{\tau} of the Borel transform
  2578. of various coefficients in the system coupling matrix.}.
  2579. As explained in the
  2580. ZOOM web page document
  2581. @url{http://www.lmd.jussieu.fr/@/ZOOM/doc/@/Feedback_Gain.pdf},
  2582. this allows for the calculation of the
  2583. dynamic gain and factor of any feedback that goes through a unique
  2584. transfer variable. An example of the conclusions that can be drawn from such
  2585. an analysis is provided in the same document.
  2586. @ignore
  2587. The Borel sweep allows the numerical study of the dependency in @math{\tau}
  2588. of the Borel transform of various
  2589. coefficients in the system coupling matrix. For example, the coefficient
  2590. @math{g(\tau)} of the reduced scalar coupling matrix
  2591. associated with a transfer defines a feedback gain.
  2592. @end ignore
  2593. For linear systems -- whose GTLS are autonomous along the whole trajectory --
  2594. the @math{\tau} function of the
  2595. feedback gain is independent of the position on the system trajectory.
  2596. But in general it is dependant, and one can analyse the function
  2597. @math{g(\tau;t)} defined on a segment @math{t} of the trajectory.
  2598. The document introducing the TEF-ZOOM technique explains how a Crank-Nicolson
  2599. scheme for the time discretisation
  2600. symbolically gives the solution of the Borel transform of the system. One can
  2601. identify the @code{dt} variable with the Borel @math{\tau} within a
  2602. factor @math{2}. Hence, to numerically study the @math{\tau} dependency of
  2603. the transform of various coefficients in the system coupling matrix at one
  2604. point in time, one can calculate the Borel transform of the TLS solutions
  2605. by making a time-step sweep.
  2606. The function @math{g(\tau;t)} is simply output for the feedback gain
  2607. attached to a unique @code{ff(k)} transfer variable.
  2608. All the relevant informations should be entered in the @file{zinit} sequence.
  2609. @menu
  2610. * Specifying the Borel sweep::
  2611. * Borel sweep results::
  2612. @end menu
  2613. @node Specifying the Borel sweep
  2614. @subsection Specifying the Borel sweep
  2615. @vindex ZBorel
  2616. First of all the logical flag @code{ZBorel} should be raised:
  2617. @example
  2618. ZBorel=.true.;
  2619. @end example
  2620. @vindex index_ff_gain
  2621. The index of the studied transfer is given in the @code{index_ff_gain}
  2622. variable
  2623. @example
  2624. index_ff_gain=7;
  2625. @end example
  2626. At each time step a Borel sweep may be performed. The time steps of interest
  2627. are
  2628. specified with three variables, one for the first step, one for the last step
  2629. and one for the number of steps between two Borel sweeps:
  2630. @vtable @code
  2631. @item istep_B_deb
  2632. First time step for the Borel sweep.
  2633. @item istep_B_fin
  2634. Last time step for the Borel sweep.
  2635. @item istep_B_inc
  2636. Number of time steps between Borel sweeps.
  2637. @end vtable
  2638. In the following examples Borel sweeps are performed from the
  2639. time step 1000 up to the time step 1200, with a sweep at each time step:
  2640. @example
  2641. istep_B_deb=1000;
  2642. istep_B_fin=1200;
  2643. istep_B_inc=1;
  2644. @end example
  2645. For each Borel sweep, the range of the @math{\tau} variable should be
  2646. set. As this is a multiplicative variable the initial value, a multiplicative
  2647. factor and the number of values are to be given.
  2648. @vtable @code
  2649. @item tau_B_ini
  2650. Initial value for @math{\tau}.
  2651. @item tau_B_mult
  2652. Multiplicative factor for sweep in @math{tau}.
  2653. @item itau_max
  2654. Number of @math{\tau} values.
  2655. @end vtable
  2656. For example, in the following, at each time step, the Borel
  2657. transform will be computed for @math{\tau} values
  2658. starting at @math{0.2} and then multiplied a hundred times by @math{\sqrt{\sqrt{2}}}
  2659. @example
  2660. tau_B_ini=0.2;
  2661. tau_B_mult=sqrt(sqrt(2.));
  2662. itau_max=100;
  2663. @end example
  2664. When the initial value of @math{\tau} is set to a negative value
  2665. (@i{i.e.} @code{tau_B_ini=-0.2;}),
  2666. the Borel sweep will first be applied with @code{itau_max} negative values
  2667. for @code{-0.2}, @code{tau_B_mult*(-0.2)},..., then for the zero value,
  2668. and finally for the symetric positive values, resulting in @code{2*itau_max+1}
  2669. values for @math{\tau}.
  2670. The whole example reads
  2671. @example
  2672. ! -------------------
  2673. ! Feedback gain
  2674. ! Borel
  2675. ! -------------------
  2676. ZBorel=.true.;
  2677. if ZBorel
  2678. < istep_B_deb=1000;
  2679. istep_B_fin=1200;
  2680. istep_B_inc=1;
  2681. ;
  2682. index_ff_gain=7;
  2683. tau_B_ini=0.2;
  2684. tau_B_mult=sqrt(sqrt(2.));
  2685. itau_max=100;
  2686. z_pr/Borel/:tau_B_mult,tau_B_ini*(tau_B_mult)**itau_max;
  2687. >;
  2688. @end example
  2689. @findex zborel for
  2690. Instead of using the index of the transfer in @code{index_ff_gain} it is
  2691. possible to specify the name of the transfer.@c , whenever
  2692. @c the symbolic model description is used (@pxref{Symbolic model description}).
  2693. In that case the transfer is specified
  2694. by the @code{zborel for} macro. For example if the transfer selected for the
  2695. feedback gain computation is @var{b_transfer}, it can be selected
  2696. with:
  2697. @example
  2698. zborel for: @var{b_transfer};
  2699. @end example
  2700. @node Borel sweep results
  2701. @subsection Borel sweep results
  2702. @cindex Borel sweep results
  2703. @cindex results, Borel sweep
  2704. @cindex Borel sweep graphics
  2705. @cindex graphics, Borel sweep
  2706. The file @file{tau_Borel.data} gives the @math{\tau} values of the @var{tau} sweep,
  2707. and the file @file{gains.data} records the feedback gain function values of
  2708. @math{g(\tau)}, with
  2709. one line for each sweep along the trajectory. In the 1.01 version, a new
  2710. feature is also provided giving the poles and residuals of the Borel
  2711. transform in the file @file{vpgains.data}. Consult the subroutine
  2712. @code{Boreleig}
  2713. for (not definitive) output description.
  2714. One can easily obtain the surface contours of @math{g(t,\tau)} using
  2715. the Fortran program provided as @file{gains.f} and its compilation shell
  2716. @file{gains.xqt},
  2717. that builds 2D histograms for PAW, in which one uses the
  2718. @file{borels.kumac} provided kumac.
  2719. @node Stability of fastest modes
  2720. @section Stability analysis of fastest modes
  2721. @cindex SVD
  2722. @cindex Singular Value Decomposition
  2723. @cindex state matrix
  2724. @cindex @file{sltc.exe}
  2725. The preceding analyses are done along with a simulation. One has also the
  2726. possibility of using in a more classical fashion the state advance matrix
  2727. @math{A_{st}}, after the end of the simulation. Code to perform the
  2728. @acronym{SVD, Singular Value Decomposition} of the state matrix @math{A_{st}}
  2729. and also of @math{A_{st} + A_{st}^\dagger} is provided with @Minik{}.
  2730. The singular elements of these two matrices correspond to the most
  2731. rapid modes of instability of the perturbed system.
  2732. The Singular value decomposition of a matrix is noted
  2733. @tex
  2734. $$
  2735. U w V^\dagger
  2736. $$
  2737. @end tex
  2738. @ifnottex
  2739. @noindent @math{U w V^t}
  2740. @end ifnottex
  2741. An executable file, @file{sltc.exe} is generated and running this file will
  2742. produce the corresponding results.
  2743. @menu
  2744. * SVD with cmz::
  2745. * SVD with make::
  2746. * SVD run and output::
  2747. @end menu
  2748. @node SVD with cmz
  2749. @subsection Singular Value Decomposition with cmz
  2750. @cindex @command{smod}
  2751. The cmz macro @code{smod SLTC} prepares a main program
  2752. (@file{circul} of +PATCH SLTC), provided as a base for user's own analysis,
  2753. in the directory @file{sltc/}.
  2754. @node SVD with make
  2755. @subsection Singular Value Decomposition with make
  2756. @cindex @file{Makefile.sltc}
  2757. To compile the singular value decomposition executable with @command{make} you
  2758. can do
  2759. @example
  2760. make sltc.exe
  2761. @end example
  2762. If you want to have a separate directory for the SVD, you should copy
  2763. the sequence @file{dimetaphi.inc} (or make a link to that file) to the
  2764. directory. You should also copy the file @file{Makefile.sltc} from the
  2765. @file{template/} directory in this directory, rename it @file{Makefile}
  2766. and set the @Minik{} directory path in the
  2767. @code{miniker_dir} variable. For
  2768. example, if the @Minik{} directory is in @file{/u/src/mini_ker}:
  2769. @example
  2770. miniker_dir = /u/src/mini_ker
  2771. @end example
  2772. @node SVD run and output
  2773. @subsection Singular Value Decomposition run and output
  2774. @cindex SVD run
  2775. @cindex run, SVD
  2776. @cindex SVD output
  2777. @cindex output, SVD
  2778. @cindex @file{sltc.exe}
  2779. @cindex @file{title.tex}, SVD
  2780. @cindex @file{aspha.data}, SVD
  2781. As it is, the @file{sltc.exe} executable generated by the compilation
  2782. determines the SVD. This program requires @file{title.tex} (@pxref{Title file}) to
  2783. transmit a title for output and graphics, and @file{aspha.data}
  2784. (@pxref{Simulation and output,,Running a simulation and using the output})
  2785. to access the
  2786. state matrix. To get access to these files (in case they are not in the current
  2787. directory) it is possible to make a link to
  2788. the corresponding files in the model directory. Once it is done
  2789. the program may be run:
  2790. @example
  2791. ./sltc.exe
  2792. @end example
  2793. The files @file{u.data}, @file{w.data}, and @file{v.data} holds the singular elements
  2794. for @math{A_{st}} (@math{U}, @math{w} and @math{V}),
  2795. and @file{us.data}, @file{ws.data}, and @file{vs.data}
  2796. holds the singular elements of @math{A_{st} + A_{st}^\dagger}.
  2797. The corresponding macros @samp{.kumac} for PAW@footnote{Explanation in
  2798. the research paper about SLTC (Al1 2003) available on request.}
  2799. are also generated.
  2800. @node Generalized TLS
  2801. @section Generalized linear tangent system analysis
  2802. @cindex Generalized linear tangent system
  2803. @cindex GTLS
  2804. @cindex propagator
  2805. @cindex Lyapunov exponents
  2806. @cindex @file{sltcirc.exe}
  2807. The state matrix @math{A_{st}} may also be used to compute the
  2808. GTLS propagator (or state transition matrix applied to perturbation), after the simulation.
  2809. The algorithm is a finite product of
  2810. 5th order development of
  2811. @math{\Phi(t+\delta t,t)=\exp{A_{st} \delta t}}.
  2812. Numerous element of analysis are given, in particular the determination
  2813. of the Lyapunov exponents of the system.
  2814. An executable file, @file{sltcirc.exe} is generated and running this file will
  2815. produce the corresponding results.
  2816. @menu
  2817. * GTLS with cmz::
  2818. * GTLS with make::
  2819. * GTLS run and output::
  2820. @end menu
  2821. @node GTLS with cmz
  2822. @subsection Generalized tangent linear system with cmz
  2823. @cindex @command{smod}
  2824. The cmz macro @code{smod SLTCIRC} prepares a main program
  2825. (@file{circule} of +PATCH SLTCIRC), in the directory @file{sltcirc/}.
  2826. @node GTLS with make
  2827. @subsection Generalized tangent linear system with make
  2828. @cindex @file{Makefile.sltcirc}
  2829. To compile the GTLS analysis executable with @command{make} you
  2830. can do
  2831. @example
  2832. make sltcirc.exe
  2833. @end example
  2834. If you want to have a separate directory for the GTLS analysis, you should copy
  2835. the sequence @file{dimetaphi.inc} (or make a link to that file) to the
  2836. directory. You should also copy the file @file{Makefile.sltcirc} from the
  2837. @file{template/} directory in this directory and rename it @file{Makefile}
  2838. and set the @Minik{} directory path in the @code{miniker_dir} variable.
  2839. @node GTLS run and output
  2840. @subsection Generalized tangent linear system analysis run and output
  2841. @cindex GTLS run
  2842. @cindex run, GTLS
  2843. @cindex GTLS output
  2844. @cindex output, GTLS
  2845. @cindex @file{sltcirc.exe}
  2846. @cindex @file{title.tex}, GTLS
  2847. @cindex @file{dres.data}, GTLS
  2848. @cindex @file{aspha.data}, GTLS
  2849. The @file{sltcirc.exe} executable generated by the compilation
  2850. computes the elements of analysis of the system. This program requires
  2851. @file{title.tex} to
  2852. transmit a title for output and graphics (@pxref{Title file}),
  2853. @file{aspha.data} to access the
  2854. state matrix and @file{dres.data}, because time-step can be changed along the
  2855. simulation
  2856. (@pxref{Simulation and output,,Running a simulation and using the output})
  2857. @footnote{cf our research texts about propagator analyses in
  2858. SLTC, and ``les Gains sur champs (Al1 2003-2004)''}. To get access to these files
  2859. (in case they are not in the current
  2860. directory) it is possible to make a link to
  2861. the corresponding files in the model directory. Once it is done
  2862. the program may be run:
  2863. @example
  2864. ./sltcirc.exe
  2865. @end example
  2866. The following table gives the correspondence between variable name,
  2867. result file and ntuple number, with a short explanation:
  2868. @multitable {@code{lwr(.,.)}} {@file{wlphit.data}} {ntuple} {eigen factors of @math{w} in the SVD of @math{\Phi}}
  2869. @headitem var @tab file @tab ntuple @tab explanation
  2870. @item @code{p(.,.)} @tab @file{phit.data} @tab 55 @tab propagator from 0 to @math{t}, @math{\Phi(t,0)}
  2871. @item @code{up(.,.)} @tab @file{uphit.data} @tab 50 @tab Left singular vectors @math{U} in the SVD of @math{\Phi}
  2872. @item @code{wp(.)} @tab @file{wphit.data} @tab 51 @tab singulat values @math{w} in the SVD of @math{\Phi}
  2873. @item @code{vp(.,.)} @tab @file{vphit.data} @tab 52 @tab Right Singular Vectors @math{V} in the SVD of @math{\Phi}
  2874. @item @code{wr(.)} @tab @file{wr.data} @tab 53 @tab real part of eigen values of @math{\Phi(t,0)}
  2875. @item @code{wi(.)} @tab @file{wi.data} @tab 54 @tab imaginary part of eigen values of @math{\Phi(t,0)}
  2876. @item @code{lwp(.)} @tab @file{lwphit.data} @tab 67 @tab Lyapunov exponents
  2877. @c @item @code{lwr(.,.)} @tab @file{lwr.data} @tab 68 @tab
  2878. @c @item @code{lwi(.,.)} @tab @file{lwi.data} @tab 69 @tab
  2879. @c @item @code{} @tab @file{.data} @tab @tab
  2880. @end multitable
  2881. @ignore
  2882. ntuple # var name
  2883. ---------------------------------------------------------------------
  2884. uphit.data 50 up SVD of Phit
  2885. wphit.data 51 wp = U[w_diag]V'
  2886. vphit.data 52 vp
  2887. wr.data 53 wr[i] eigen factors
  2888. wi.data 54 wi real and imaginary parts
  2889. ---------------------------------------------------------------------
  2890. lphit.data 65 lp[ij] 1/t Ln Phit + 1/LnDetPhit,<TrA>
  2891. ulphit.data 60 ulp SVD of 1/t LnPhit
  2892. wlphit.data 61 wlp
  2893. vlphit.data 62 vlp
  2894. wrl.data 63 wrl[i] VP of 1/t LnPhit
  2895. wil.data 64 wil
  2896. ---------------------------------------------------------------------
  2897. lwphit.data 67 lwp[i] 1/t Ln w : Phit=UwV'
  2898. lwr.data 68 lwr 1/t VP de UV'
  2899. lwi.data 69 lwi
  2900. ---------------------------------------------------------------------
  2901. The bloc 60-65 contains the elements of analysis of the propagator when
  2902. it is assimilated to the autonomous system such that @math{\Phi(t,0)=\exp{At}},
  2903. where @math{A} is a state matrix of an autonomous system giving the same
  2904. sensitivity to perturbation at time @math{t}.
  2905. Last bloc 67-69 gives the Lyapunov exponents along with their
  2906. corresponding eigen vectors.
  2907. @end ignore
  2908. @node Advanced use of @Minik{} with make
  2909. @chapter Advanced use of @Minik{} with make
  2910. @menu
  2911. * Make variables::
  2912. * Rules::
  2913. * Linking rule::
  2914. @end menu
  2915. @node Make variables
  2916. @section Make variables
  2917. @cindex @file{Makefile.miniker}
  2918. The @file{Makefile.miniker} Makefile provided in the
  2919. distribution should be included as it defines a lot of important
  2920. variables and rules.
  2921. The following make variables can be set by the user:
  2922. @table @code
  2923. @item miniker_dir
  2924. that variable should hold the @Minik{} sources directory. If you installed
  2925. @Minik{} that variable should be set to @file{$(includedir)/mini_ker}.
  2926. If you use the sources right from the sources directory it should be set to
  2927. the sources package directory.
  2928. @item MTNDIRS
  2929. This variable can hold a @samp{:} delimited list of directories that will
  2930. be searched for mortran include files.
  2931. @item CMFDIRS
  2932. This variable can hold a @samp{:} delimited list of directories that will
  2933. be searched for cmz directive include files.
  2934. @item SEL
  2935. This variable holds a @samp{,} delimited list of select flags, for example
  2936. @code{monitor}, @code{grid1d}, @code{debug}.
  2937. @item LDADD
  2938. This variable can be used to add libraries flags and files. It is used in
  2939. the default linking command/rule.
  2940. @item miniker_user_objects
  2941. This variable should hold a space separated list of additional object files
  2942. to be linked with the model and helper object files.
  2943. @item CAR2TXTFLAGS
  2944. cmz directives preprocessor flag.
  2945. @item kalman
  2946. This variable should be set to 1 if you want to use the kalman filter
  2947. (@pxref{Kalman filter}).
  2948. @item double
  2949. This variable should be set to 1 if you want to have a double precision
  2950. code (@pxref{Double precision}).
  2951. @end table
  2952. The following variables are allready set and may be used
  2953. (some are set by ./configure see @ref{Configuration}):
  2954. @table @code
  2955. @item miniker_principal_objects
  2956. The list of object files needed for the model build, together with some
  2957. helper object files often used but not strictly required for the linking.
  2958. @item DEPDIR
  2959. The name of a hidden directory containing the dependencies computed
  2960. for the main mortran files.
  2961. @itemx F77
  2962. @itemx FC
  2963. @itemx FFLAGS
  2964. @item LDFLAGS
  2965. Compiler and linker related variables set by ./configure.
  2966. @item LIBS
  2967. This variable should hold the link flags and files required to build
  2968. @Minik{}, set by ./configure.
  2969. @item CAR2TXT
  2970. @itemx MORTRAN
  2971. @itemx MTNFLAGS
  2972. @itemx MTNDEPEND
  2973. Preprocessor and preprocessor flags, set by ./configure.
  2974. @end table
  2975. @node Rules
  2976. @section Rules
  2977. The following rules are defined in the @file{Makefile.miniker} file.
  2978. @table @code
  2979. @item miniker-clean
  2980. remove the fortran files generated from the mortran files. Remove
  2981. the object files.
  2982. @item miniker-mtn-clean
  2983. remove the mortran files generated from the files with cmz directives.
  2984. @item
  2985. Various rules to preprocess files with cmz directives and mortran files and
  2986. to compile fortran files.
  2987. @end table
  2988. If the user needs a mortran main file, he may take advantage of the rule
  2989. used to compute the dependencies of a mortran file. If the file is called,
  2990. say, @file{mtnfile.mtn} leading to @file{mtnfile.f}, the following include
  2991. should lead to the automatic creation, updating and inclusion of a
  2992. file describing the dependencies of @file{mtnfile.mtn} in the
  2993. @file{Makefile}:
  2994. @example
  2995. include $(DEPDIR)/mtnfile.Pf
  2996. @end example
  2997. @node Linking rule
  2998. @section Linking rule
  2999. The rule used for the linking of the model file is not in the
  3000. @file{Makefile.miniker} file but
  3001. should be provided in the user @file{Makefile} for more flexibility.
  3002. The default rule
  3003. uses the variables @code{miniker_user_objects} for additional object files
  3004. and @code{LDADD} for additionnal linking flags and files, those
  3005. variables are there to be changed by the user.
  3006. The object files required by the @Minik{} code are in the make variable
  3007. @code{miniker_principal_objects}, this variable is also used.
  3008. The value of the variables @code{FC}
  3009. for the Fortran compiler, @code{FFLAGS} for the Fortran compiler
  3010. flags and @code{LDFLAGS} for the linker flags should be set to right
  3011. values; @code{LIBS} should also be right and hold the link flags and link
  3012. files required to compile the @Minik{} model. These variables are
  3013. set by by @command{./configure} during configuration (@pxref{Configuration})
  3014. and used in the default rule:
  3015. @verbatim
  3016. $(model_file): $(miniker_user_objects) $(miniker_principal_objects)
  3017. $(FC) $(FFLAGS) $(LDFLAGS) $^ $(LDADD) $(LIBS) -o $@
  3018. @end verbatim
  3019. In case this isn't right it may be freely changed. You should certainly
  3020. refer to the @ref{Top,,Top,make,GNU Make Manual} manual to understand what
  3021. that rule exactly means and make your own.
  3022. @node Concepts index
  3023. @unnumbered Concepts index
  3024. @printindex cp
  3025. @node Variables macros and functions index
  3026. @unnumbered Variables, macros and functions index
  3027. @printindex vr
  3028. @node Installation
  3029. @appendix Installation
  3030. @menu
  3031. * Programming environments::
  3032. * Common requisites::
  3033. * @Minik{} with cmz::
  3034. * @Minik{} with make::
  3035. @end menu
  3036. @node Programming environments
  3037. @appendixsec Programming environments
  3038. @cindex Programming environments
  3039. @Minik{} is not a traditionnal software in that it isn't a library
  3040. or an interpreter but rather a set of source and macro file that
  3041. combines with the user model code and enable
  3042. to build a binary program corresponding with the model. It
  3043. requires a build environment with a preprocessor, a compiler
  3044. and facilities that automate these steps.
  3045. Two different environment are proposed. One use
  3046. @command{cmz} (@url{http://wwwcmz.web.cern.ch/@/wwwcmz/index.html}),
  3047. while the other is based on @command{make}. Other libraries
  3048. are needed, the CERN Program Library (cernlib) and lapack.
  3049. @node Common requisites
  3050. @appendixsec Common requisites
  3051. @cindex cernlib
  3052. @cindex lapack
  3053. Whatever method is used a fortran 77 compiler is required. The compilers
  3054. that have been used so far are g77, gfortran and the sun solaris compiler.
  3055. When usng CMZ, the CERN Program Library, available at
  3056. @url{http://wwwasd.web.cern.ch/@/wwwasd/cernlib/}, has to be installed.
  3057. With make, internal source files copied from the cernlib may be used instead
  3058. but then some examples won't be available, since they rely on some
  3059. mathematical functions provided by the CERN library.
  3060. On windows, in case you want to use the compiler from the GNU compiler
  3061. collection with cygwin or MINGW/MSYS you can use the binaries provided at
  3062. @url{http://zyao.home.cern.ch/@/zyao/cernlib.html}.
  3063. On Mac OS X, the cernlib provided by fink (package @code{cernlib-devel})
  3064. can be used.
  3065. You should also have LAPACK, available at @url{http://www.netlib.org/@/lapack/}.
  3066. LAPACK can also be installed as part of the CERN Library or as part of
  3067. the @uref{ATLAS,http://math-atlas.sourceforge.net/} implementation.
  3068. On most linux distributions a lapack package is available.
  3069. On Mac OS X, the ATLAS implementation provided by fink or the frameworks
  3070. from Xcode can be used.
  3071. @node @Minik{} with cmz
  3072. @appendixsec @Minik{} with cmz
  3073. @cindex @file{mini_ker.cmz}
  3074. @cindex @file{selseq.kumac}
  3075. First of all you have to get the cmz file @file{mini_ker.cmz} and put it
  3076. in a directory. In that same directory you should create a directory for
  3077. each of your models. In the model directory you should copy the file
  3078. @file{selseq.kumac} available with @Minik{}, and create your own cmz
  3079. file for your model, called for example @file{mymodel.cmz}. You should also
  3080. have installed the kumac macro files
  3081. handling mortan compilation, the associated shell scripts and the mortran
  3082. preprocessor.
  3083. @node @Minik{} with make
  3084. @appendixsec @Minik{} with make
  3085. @menu
  3086. * Additional requirements::
  3087. * Configuration::
  3088. * Installation with make::
  3089. @end menu
  3090. @node Additional requirements
  3091. @appendixsubsec Additional requirements for @Minik{} with make
  3092. @cindex @command{mortran}, with make
  3093. @cindex requirements, with make
  3094. The package has been tested with GNU @command{make} and solaris
  3095. @command{make}.
  3096. Suitable preprocessors should also be installed. Two preprocessors are
  3097. required, one that preprocess the cmz directives, and a mortran
  3098. preprocessor. A cmz directives processor written in @command{perl},
  3099. is distributed in the @command{car2txt} package available at
  3100. @value{myurl}. A @command{mortran}
  3101. package with a command able to preprocess a mortran file given on
  3102. the command line with a syntax similar with the @command{cpp} command line
  3103. syntax is also required.
  3104. Such a mortran is available at @value{myurl}.
  3105. @c All the @command{make} commands are not suitable, for example the distribution
  3106. @c doesn't work with solaris @command{make}. GNU @command{make} works, however,
  3107. @c and should be available on most platforms, often called @command{gmake}.
  3108. @node Configuration
  3109. @appendixsubsec Configuration
  3110. @cindex configuration of source
  3111. The package is available at @value{myurl}. It is
  3112. available as a compresssed tar archive. On UNIX, with GNU @command{tar} it
  3113. may be unpacked using
  3114. @example
  3115. $ tar xzvf mini_ker-@value{VERSION}.tar.gz
  3116. @end example
  3117. The detection of the compiler, the preprocessors (car2txt and mortran),
  3118. and the libraries are performed by the configure script. This script
  3119. sets the
  3120. apropriate variables in makefiles. It can be run with:
  3121. @example
  3122. $ cd mini_ker-@value{VERSION}
  3123. $ ./configure
  3124. @end example
  3125. If the output of @command{./configure} doesn't show any error it means that
  3126. all the components are here. It is possible to give @command{./configure}
  3127. switches and also specify environment variables (see also
  3128. @command{./configure --help}):
  3129. @table @code
  3130. @item --disable-cernlib
  3131. Use the internal cernlib source files, even if a cernlib is detected.
  3132. @item --with-static-cernlib
  3133. This command line switch forces a static linking with the cernlib (or a dynamic linking
  3134. if set to no).
  3135. @item --with-cernlib
  3136. This command line switch can be used to specify the cernlib location
  3137. (if not detected or you want to use a specific cernlib).
  3138. @item --with-blas
  3139. @itemx --with-lapack
  3140. With this command switch, you can specify the location of the blas and lapack
  3141. libraries.
  3142. For example, on mac OS X this can be used to specify the blas and lapack from
  3143. the Apple frameworks:
  3144. @example
  3145. ./configure \
  3146. --with-blas=/System/Library/Frameworks/vecLib.framework/versions/A/vecLib \
  3147. --with-lapack=/System/Library/Frameworks/vecLib.framework/versions/A/vecLib
  3148. @end example
  3149. @item F77
  3150. @itemx FC
  3151. @itemx FFLAGS
  3152. @itemx LDFLAGS
  3153. Classical compiler, compiler flags and linker flags.
  3154. @item MORTRAN
  3155. This environment variable holds the mortran preprocessor command
  3156. (default is @command{mortran}).
  3157. @item MTNFLAGS
  3158. This environment variable holds command line arguments for the mortran
  3159. preprocessor. It is empty in the default case.
  3160. @item MTN
  3161. This environment variable may be used to specify the mortran executable
  3162. name and/or path, it should be used by the @command{mortran} commmand.
  3163. (default is empty, which leads to a mortran executable called @command{mtn}).
  3164. @item MTNDEPEND
  3165. This environment variable may be used to specify the mortran dependencies
  3166. checker executable. It should be used by the @command{mortran} commmand.
  3167. (default is empty, which leads to a mortran dependencies checker
  3168. called @command{mtndepend}).
  3169. @end table
  3170. After a proper configuration, if @command{make} is run then the example
  3171. models should be build. You have to perform the configuration only once.
  3172. @node Installation with make
  3173. @appendixsubsec Installation with make
  3174. @cindex installation with make
  3175. @Minik{} can be installed by running
  3176. @example
  3177. make install
  3178. @end example
  3179. It should copy the sources
  3180. and the @file{Makefile.miniker} file in
  3181. a @file{mini_ker} directory in the @code{$(includedir)} directory, and
  3182. copy the templates in @file{$(datadir)/mini_ker}. The default for
  3183. @code{$(includedir)} is @file{/usr/local/include} and the default for
  3184. @code{$(datadir)} is @file{/usr/local/share}, these defaults may be
  3185. changed by @command{./configure} switches @samp{--prefix},
  3186. @samp{--includedir} and @samp{--datadir}. See @command{./configure --help}
  3187. and the @file{INSTALL} file for more informations. The helper script
  3188. @file{start_miniker} should also be installed.
  3189. The installation is not required to use comfortably @Minik{}. Indeed
  3190. the only thing that changes with the sources and the @file{Makefile.miniker}
  3191. directory location is the @code{miniker_dir} variable in a
  3192. project @code{Makefile}.
  3193. @node Cmz directives reference
  3194. @appendix Cmz directives reference
  3195. The cmz directives are described together with the other
  3196. features of cmz in the cmz manual at
  3197. @url{http://wwwcmz.web.cern.ch/wwwcmz/}, the important ones are
  3198. nevertheless recalled here,
  3199. especially for those that use make and don't need the whole
  3200. features of cmz.
  3201. After the description of the generic features, we turn
  3202. to the cmz directive of interest.
  3203. There are three kinds of cmz directives that are of use
  3204. within @Minik{}: one kind
  3205. that introduce files, the other for conditionnal compilation and
  3206. the third for sequence inclusion.
  3207. @menu
  3208. * Cmz directives general syntax::
  3209. * Conditional expressions::
  3210. * File introduction directives::
  3211. * Conditional directives::
  3212. * File inclusion directive::
  3213. * The self directive::
  3214. @end menu
  3215. @node Cmz directives general syntax
  3216. @appendixsec Cmz directives general syntax
  3217. The cmz directives always begin with a @samp{+} in the first column,
  3218. optionnaly followed by any number of @samp{_} that may be used for
  3219. indentation, then the directive label, case insensitive, followed
  3220. by the directive arguments separated by @samp{,}. The arguments
  3221. are also case insensitive.
  3222. Optional spaces may be around directive arguments.
  3223. An optionnal @samp{.} ends the directive
  3224. arguments and begin a comment, everything that follows that @samp{.} is
  3225. ignored.
  3226. @node Conditional expressions
  3227. @appendixsec Conditional expressions
  3228. A directive argument common to all the directives is the conditionnal
  3229. expression. A conditionnal expression may be true or false, it is a
  3230. combination of select flags. the select flags are combined with
  3231. logical operators. A
  3232. select flag itself is true if it was selected. A select flag @var{selflag}
  3233. is selected by using the @code{sel @var{selflag}} instruction in cmz. It is
  3234. selected by passing the @code{-D @var{selflag}} command line switch to
  3235. the call of the cmz directives preprocessor when using make.
  3236. A @samp{-} negates
  3237. the expression that follows. Parenthesis @samp{(} and @samp{)} are used
  3238. for the grouping of subexpressions. @samp{|} and @samp{,} are for the
  3239. boolean or: an expression with a or is true
  3240. if the expression on the left or the expression on the right of the or
  3241. is true.
  3242. @samp{&} is for the boolean and: an expression with an and is true if
  3243. the expression on the left and the expression on the right are true.
  3244. The grouping is left to right when there is no parenthesis, with or and
  3245. @samp{&} having the same precedence. Therefore
  3246. @example
  3247. a&b|c @equiv{} (a&b)|c
  3248. a|b&c @equiv{} (a|b)&c
  3249. a|b&c is not a|(b&c)
  3250. a&b|c is not a&(b|c)
  3251. @end example
  3252. @node File introduction directives
  3253. @appendixsec File introduction directives
  3254. A file (or sequence) introduction directive appears at the beginning
  3255. of the file. There are two different directives, one is @code{DECK}
  3256. for normal files, the other is @code{KEEP} for include files (sequences).
  3257. The first argument is the name of the file. The file name may not be larger
  3258. than 32 characters and is converted to lower case in the general case.
  3259. The optionnal following arguments may be
  3260. of 2 type (and may be mixed, separated by @samp{,}):
  3261. @table @asis
  3262. @item conditional
  3263. A conditionnal is introduced by @code{IF=} followed by a conditionnal
  3264. expression described in
  3265. @ref{Conditional expressions}. The
  3266. file is preprocessed if the conditionnal expression is true.
  3267. @item language specification
  3268. A language specification is introduced by a @code{T=}. The most
  3269. common languages are @samp{mtn} for the mortran, @samp{ftn} for
  3270. fortran not preprocessed, @samp{f77} for preprocessed fortran,
  3271. @samp{c} for the c language and @samp{txt} for text files.
  3272. In general the language of the file determines the name of files
  3273. the preprocessed file is extracted to, the comment style and
  3274. the command for inclusions.
  3275. @end table
  3276. It is a common practice to have wrong language type in @code{KEEP}
  3277. as the language may be determined from the @code{DECK} that include
  3278. them with cmz, or from their file name with make. This is not recommended
  3279. and considered a bad practice.
  3280. Such a directive will always appear in cmz, as it is built-in. It
  3281. is recommended to have one when using make too, even though it is not
  3282. required in most cases. Indeed make uses the file name directly
  3283. and finds the language and file type by looking at the file extension.
  3284. make should then pass the language type with a
  3285. @code{--lang @var{lang}} command
  3286. line switch when calling the cmz directives preprocessor.
  3287. With make, the convention is to have @samp{cm} added before the normal
  3288. file suffix and after the @samp{.}. The table @ref{tab:cmfile_suffix}
  3289. shows the matching between suffixes, file type and file language.
  3290. For example, a file beginning with
  3291. @verbatim
  3292. +Deck, subroutine_foo, If=monitor&-simple, T=f77.
  3293. @end verbatim
  3294. is a main preprocessed fortran file that will only be generated if
  3295. @samp{monitor} is selected and @samp{simple} is not selected. The
  3296. file to be preprocessed by make should have the @samp{.cmF} suffix,
  3297. and be called @file{subroutine_foo.cmF}.
  3298. A file beginning with
  3299. @verbatim
  3300. +KEEP,inc_common,If=monitor|interface,T=mtn
  3301. @end verbatim
  3302. is an mortran include file that should be processed only if @samp{monitor}
  3303. or @samp{interface} is selected. The file to be preprocessed by make
  3304. should have the @samp{cmmti} suffix and be called @file{inc_common.cmmti}.
  3305. The resulting file when make is used will be called @file{inc_common.mti}.
  3306. @node Conditional directives
  3307. @appendixsec Conditional directives
  3308. Conditional directives may be used to conditionnaly skip blocks of
  3309. code. There are 4 conditional directives: @code{if}, @code{elseif},
  3310. @code{else} and @code{endif}. @code{+if} begins a conditional directives
  3311. sequence, with argument a conditional expression. If the expression is
  3312. true the block of code following the @code{+if} is output in the
  3313. resulting file, up to another conditional directive, if it is false
  3314. the code block is skipped. If the
  3315. expression is false and the following conditional directive is
  3316. @code{+elseif}, the same procedure is followed with the argument of
  3317. @code{+elseif}
  3318. which is also a conditionnal expression. More than one @code{+elseif}
  3319. may follow a @code{+if}. If a @code{+if} or @code{+elseif} expression
  3320. is true the following
  3321. code block is output and all
  3322. the following @code{+elseif} code blocks are skipped. If all the @code{+if}
  3323. and @code{+elseif} expressions are false and
  3324. the following coditionnal
  3325. directive is @code{+else} then the block following the
  3326. @code{+else} is output. If a previous expression was true the
  3327. code block following the @code{+else} is skipped. The last code block
  3328. is closed by @code{+endif}.
  3329. Conditionnal directives may be nested, a @code{+if} begins a deeper
  3330. conditionnal sequences directives that is ended by the corresponding
  3331. @code{+endif}.
  3332. The simplest example is:
  3333. @verbatim
  3334. some code;
  3335. +IF,monitor
  3336. code output only if monitor is true;
  3337. +ENDIF
  3338. @end verbatim
  3339. If @samp{monitor} is selected, the @code{+if} block is output, it leads to
  3340. @verbatim
  3341. some code;
  3342. code output only if monitor is true;
  3343. @end verbatim
  3344. If @samp{monitor} isn't selected the @code{+if} block is skipped, it leads to
  3345. @verbatim
  3346. some code;
  3347. @end verbatim
  3348. An example with @code{+else} may be:
  3349. @verbatim
  3350. +IF,double
  3351. call dmysub(eta);
  3352. +ELSE
  3353. call smysub(eta);
  3354. +ENDIF
  3355. @end verbatim
  3356. If @samp{double} is selected the code output is @code{call dmysub(eta);},
  3357. if @samp{double} isn't selected the code output is @code{call dmysub(eta);}.
  3358. Here is a self explanatory example of use of @code{+elseif}:
  3359. @verbatim
  3360. +IF,monitor
  3361. code used if monitor is selected;
  3362. +ELSEIF,kalman
  3363. code used if kalman is selected and monitor is not;
  3364. +ELSE
  3365. code used if kalman and monitor are not selected;
  3366. +ENDIF
  3367. @end verbatim
  3368. And last an example of nested conditional directives:
  3369. @verbatim
  3370. +IF,monitor
  3371. code used if monitor is selected;
  3372. +_IF,kalman. deep if
  3373. code used if monitor and kalman are selected;
  3374. +_ELSE. deep else
  3375. code used if monitor is selected and kalman is not;
  3376. +_ENDIF. end the deep conditionnals sequence
  3377. +ELSE
  3378. code used if monitor is not selected;
  3379. +_IF,kalman
  3380. code used if monitor is not selected but kalman is;
  3381. +_ELSE
  3382. code used if monitor and kalman are not selected;
  3383. +_ENDIF
  3384. other code used if monitor is not selected;
  3385. +ENDIF
  3386. @end verbatim
  3387. @node File inclusion directive
  3388. @appendixsec File inclusion directive
  3389. The file (sequence) inclusion directive is @code{seq}. The argument of
  3390. @code{seq} is an include files @samp{,} separated list. The include
  3391. files are @code{Keep} in cmz. The following optional arguments may be
  3392. mixed:
  3393. @table @asis
  3394. @item conditional
  3395. A conditionnal is introduced by @code{IF=} followed by a conditionnal
  3396. expression described in
  3397. @ref{Conditional expressions}. The
  3398. directive is ignored if the conditionnal expression is false.
  3399. @item T=noinclude
  3400. When this argument is present the text of the sequence will
  3401. always be included in the file where the @code{+seq} appears.
  3402. @end table
  3403. When there is no @code{T=noinclude} argument, the @code{+seq}
  3404. directive may be replaced with an inclusion command suitable
  3405. for the language of the file being processed, if such
  3406. command has been specified.
  3407. For example if we have the following sequence
  3408. @verbatim
  3409. +KEEP,inc,lang=C
  3410. typedef struct incstr {char* msg};
  3411. @end verbatim
  3412. And the following code in the file being processed:
  3413. @verbatim
  3414. +DECK,mainf,lang=C
  3415. +SEQ,inc
  3416. int main (int argc, char* argv) { exit(0); }
  3417. @end verbatim
  3418. the processing of @file{mainf} should lead to the file
  3419. @file{mainf.c}, containing
  3420. an include command for @file{inc}:
  3421. @verbatim
  3422. #include "inc.h"
  3423. int main (int argc, char* argv) { exit(0); }
  3424. @end verbatim
  3425. In case the @code{+seq} has the @code{T=noinclude}:
  3426. @verbatim
  3427. +DECK,mainf,lang=C
  3428. +SEQ,inc,T=noinclude
  3429. int main (int argc, char* argv) { exit(0); }
  3430. @end verbatim
  3431. The processing of @file{mainf} should lead to the file @file{mainf.c}
  3432. containing the text of @file{inc}:
  3433. @verbatim
  3434. typedef struct incstr {char* msg};
  3435. int main (int argc, char* argv) { exit(0); }
  3436. @end verbatim
  3437. @node The self directive
  3438. @appendixsec The @samp{self} directive
  3439. The @code{self} directive is an obsolete directive that may be used for
  3440. conditionnal skipping of code. For a better approach see
  3441. @ref{Conditional directives}.
  3442. The optionnal argument of @code{+SELF} is @code{If=} followed by
  3443. a conditionnal expression. If the conditionnal expression is true the
  3444. code following the directive is output, if it is false the code
  3445. is skipped up to any directive (including another @code{+SELF})
  3446. except @code{+seq}.
  3447. @ignore
  3448. @node Resolution method
  3449. @appendix Overview of resolution method
  3450. @node @Minik{} macros
  3451. @appendix @Minik{} macros
  3452. @end ignore
  3453. @node Copying This Manual
  3454. @appendix Copying This Manual
  3455. @menu
  3456. * GNU Free Documentation License:: License for copying this manual.
  3457. @end menu
  3458. @include fdl.texi
  3459. @bye