Com potser que al final lo programa funcionés???

 

Primer m’havia equivocat en un signe en el càlcul del vent geostròfic . Això feia que la sortida calculada no tingués cap sentit .

 

També vaig cometre l’error de fer servir la fòrmula  f = f0 + βy durant tot l’interval de latituds amb la f0  de 25 ºN, quan només té validesa a l’entorn d’f0  .Finalment vaig fer servir simplement : f = 2Ω sinØ.

 

Després me’n vaig adonar de que no feia bé la part del càlcul de la inversa de la laplaciana. Pel mètode de la relaxació simultània  ()  partim dels valors de la laplaciana que coneixem a cada punt ( )  per a aconseguir l’objectiu () , que no és altre que el valor de la tendència de geopotencial . Iteració per iteració , el residu () ha d’anar fent-se cada cop més petit:

 

Per a donar el número d’iteracions per bo, limitava l’acumulació dels residuos  a 0.001:

 

comptador=0;

do

{

càlcul de la inversa de la laplaciana

comptador++;

}

while ();

 

El comptador sempre donava de l’ordre de 3000 iteracions, però molt sobint acabava el procès amb nomes 150 , 200... lo que feia que el càlcul fos incorrecte. Això tenia una lògica clara: de quant en quant els residuos al sumar-se  se cancelaven , ja que tot i ser molt diferents de 0, tenien signe contrari i un valor similar. Per tant, s’havia de fer la suma dels valors absoluts del residu:

 

comptador=0;

do

{

càlcul de la inversa de la laplaciana

comptador++;

}

while ();

 

Però el procès en unes poques integracions es feia divergent. Començaven a formar-se ondulacions cada cop més grans tant entre el tàlveg i la dorsal com a la part superior esquerra del mapa.

Llavors vaig “arreglar” el programa posant un filtre “passa-baix” dels dels apunts, per a atenuar les ones curtes i així evitar les oscil.lacions. La cosa començava a anar millor. Les oscil.lacions només es donaven a la part de la baixa d’Islàndia ( a l’extrem superior esquerra del mapa) , però no a la zona entre el tàlveg i la dorsal. Hi havia encara alguna cosa que no quadrava, i havia d’estar relacionat amb el tractament dels punts de la frontera de la xarxa.

Vaig haver de programar tres tipus de derivada en  lloc de només una. Els punts que estan en els límits de la xarxa no es poden tractar com els altres. La derivada ha de ser només en una direcció: s’ha d’obtenir la diferència entre el valor del punt en qüestió i el seu veí de l’interior , però no es pot fer servir cap valor del seu veí de l’exterior ja que aquest està fora dels límits.

Lo mateix podem dir de la laplaciana. Aquesta no es pot calcular en els limits de la xarxa.

Per això vaig implementar la derivada de primer ordre als punts més exteriors, la derivada de segon ordre als veïns interiors d’aquests i la derivada de quart ordre a la resta dels punts, ja que aquesta té molta més exactitud (tal i com hem vist al problema 38).

Per a fer que els punts dels límits xarxa no queden estàtics i puguin evolucionar com els seus veïns , he aplicat una tendència de geopotencial igual a la calculada per als seus veïns interns. D’aquesta manera no hi haurà una discontinuitat en aquests punts que pugui fer que el sistema “peti”. Però el sistema seguia sense convergir. Finalment fent proves, vaig comprovar que si anulava l’efecte de curvatura , és a dir, si feia que  dx = dy  per a cada punt independentment de la latitud ( el programa tenia en compte que dx  disminuia cap al nord segons : dx = dy * cos (ø)), llavors el programa funcionava perfectament , com a mínim fins a 48 passos de 30 minuts cada un. Vist això me’n vaig adonar de que la laplaciana l’estava calculant amb el terme d*d al denominador, i això no pot ser del tot correcte, ja que la d , és a dir la distància entre dos punts de la xarxa ,no és igual en direcció N-S que en direcció E-W (no es cumpleix, com hem dit abans que dx=dy ):

 

 

Vaig provar de canviar d*d per dx*dy a veure què passava i finalment va funcionar .  No se si de casualitat o què.

De tota manera, la forma en que el procés funciona més ràpid és quan treballem amb una d constant, per tant vaig acabar “d’arreglar “ el  programa treballant només amb d=dx=dy .

 

Efecte del filtre passa-baix:

 

Aquest operador  és la peça clau del programa , ja que permet “llimar” les irregularitats a escala petita i deixar intacte el senyal a distàncies grans. Podem veure a continuació el valor de l’advecció de vorticitat absoluta en la primera integració del programa sense el filtre:

 

I amb el filtre:

 

 

 

A continuació podem veure com “funciona” el programa sense filtre, amb intervals d’integració de 15 minuts:

 

 

Però no tot lo que aporta el filtre son bondats... té l’inconvenient de que després d’aplicar-lo varies vegades, el seu efecte “llimador” se nota cada cop a distàncies més grans , de manera que el tàlveg es va reomplint i al final s’acaba transformant lo mapa en una  configuració practicament “barotròpica”( entre cometes perque jo diria que només seria barotròpica si les isobares tinguessen pendent nula ), amb les isohipses paral.leles als cercles de igual latitud ( i a les isotermes si no hi haguessin continents..):

 

 

Quin és l’orígen de la inestabilitat?

 

Com podem veure a l’anterior imatge de l’advecció de vorticitat , el filtre no pot llimar del tot les anomalies que es veuen a la part oriental de la península Ibèrica, degudes probablement al fet de que en aquells llocs és on la divergència és màxima. Per a comprovar aquesta hipòtesi, podem determinar quina és la sortida del primer pas d’integració si no negligim la divergència. És a dir si fem servir l’expressió següent:

 

 

L’equació de vorticitat quasigeostròfica anterior ens diu que la vorticitat aumentarà en un punt determinat si hi ha una advecció de vorticitat positiva, però també que disminuirà si tenim divergència positiva. Tenint en compte això , els llocs del mapa anterior on la tendència de la vorticitat està molt per damunt de la mitjana ( linies vermelles), serà perque aquesta s’ha sobreestimat al no restar el terme de la divergència , on és clarament positiva. Als llocs on la divergència sigui negativa, en canvi, aquesta tendència s’haurà infravalorat.  Aquestes variacions son amplificades i desfassades 180 º cada cop que es calcula la laplaciana:

 

 

 

Com la k = 2л / λ , si la longitud d’ona és molt curta , la k és molt gran. A mesura que el senyal inicial es va derivant, se va “embrutant” cada cop més de senyal  “d’alta freqüència “ fins que el programa se “satura”. Això passa tot i que segons el problema 38, la derivada calculada numèricament amb l’aproximació de quart ordre , es més petita respecte de la real com més petita sigui la λ respecte de Δx .

 

Si calculem la tendència de la vorticitat a cada punt tenint en compte la divergència, després del primer pas d’integració ens queda:

 

Com podem veure, no hi ha cap variació anormal a cap punt.

Si puguessim determinar d’alguna manera la tendència de la divergència a cada punt , podriem calcular la tendència del geopotencial amb més precisió, i no ens faria falta cap filtre passa-baixos.