## L. Fita, LMD October 2014. # Generation of initial conditions for an aqua-planet from the LMDZ model (r1818) 'iniacademic.F90' program # Author: Frederic Hourdin original: 15/01/93 # The forcing defined here is from Held and Suarez, 1994, Bulletin # of the American Meteorological Society, 75, 1825. from optparse import OptionParser import numpy as np from netCDF4 import Dataset as NetCDFFile import os import re import nc_var_tools as ncvar from lmdz_const import * main = 'iniaqua.py' errormsg='ERROR -- error -- ERROR -- error' warnmsg='WARNING -- warning -- WARNING -- warning' filekinds = ['CF', 'LMDZ', 'WRF'] ## e.g. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z tropo -q 2 def fxy(dx, dy): """! ! $Id: fxy.F 1403 2010-07-01 09:02:53Z fairhead $ ! c Auteur : P. Le Van c c Calcul des longitudes et des latitudes pour une fonction f(x,y) c a tangente sinusoidale et eventuellement avec zoom . c c """ fname ='fxy' #c ...... calcul des latitudes et de y' ..... #c vrlatu = np.zeros((dy+1), dtype=np.float) vyprimu = np.zeros((dy+1), dtype=np.float) vrlatu = ffy(np.arange(dy+1)*1. + 1.) vyprimu = ffy(np.arange(dy+1)*1. + 1.) vrlatv = np.zeros((dy), dtype=np.float) vrlatu1 = np.zeros((dy), dtype=np.float) vrlatu2 = np.zeros((dy), dtype=np.float) yprimv = np.zeros((dy), dtype=np.float) vyprimu1 = np.zeros((dy), dtype=np.float) vyprimu2 = np.zeros((dy), dtype=np.float) vrlatv = ffy(np.arange(dy)*1. + 1. + 0.5) vrlatu1 = ffy(np.arange(dy)*1. + 1. + 0.25) vrlatu2 = ffy(np.arange(dy)*1. + 1. + 0.75) vyprimv = fyprim(np.arange(dy)*1. + 1. + 0.5) vyprimu1 = fyprim(np.arange(dy)*1. + 1. + 0.25) vyprimu2 = fyprim(np.arange(dy)*1. + 1. + 0.75) #c #c ..... calcul des longitudes et de x' ..... #c vrlonv = np.zeros((dx+1), dtype=np.float) vrlonu = np.zeros((dx+1), dtype=np.float) vrlonm025 = np.zeros((dx+1), dtype=np.float) vrlonp025 = np.zeros((dx+1), dtype=np.float) vxprimv = np.zeros((dx+1), dtype=np.float) vxprimu = np.zeros((dx+1), dtype=np.float) vxprimm025 = np.zeros((dx+1), dtype=np.float) vxprimo025 = np.zeros((dx+1), dtype=np.float) vrlonv = ffy(np.arange(dx+1)*1. + 1.) vrlonu = ffy(np.arange(dx+1)*1. + 1. + 0.5) vrlonm025 = ffy(np.arange(dx+1)*1. + 1. - 0.25) vrlonp025 = ffy(np.arange(dx+1)*1. + 1. + 0.25) vxprimv = fxprim(np.arange(dx+1)*1. + 1.) vxprimu = fxprim(np.arange(dx+1)*1. + 1. + 0.5) vxprimm025 = fxprim(np.arange(dx+1)*1. + 1. - 0.25) vxprimp025 = fxprim(np.arange(dx+1)*1. + 1. + 0.25) return vrlatu, vyprimu, vrlatv, vyprimv, vrlatu1, vyprimu1, vrlatu2, vyprimu2, \ vrlonu, vxprimu, vrlonv, vxprimv, vrlonm025, vxprimm025, vrlonp025, vxprimp025 # From grid/fxy_sin.h # ................................................................ # ................ Fonctions in line ........................... # ................................................................ # def fy(rj,dy): """ """ val = np.arcsin(1.+2.*((1.-rj)/np.float(dy))) return val def fyprim(rj,dy): """ """ val = 1./np.sqrt((rj-1.)*(dy+1.-rj)) return val def fx(ri,dx): """ """ val = 2.*np.pi/np.float(dx) * ( ri - 0.5*np.float(dx) - 1. ) return val def fxprim(ri,dx): """ """ val = 2.*np.pi/np.float(dx) return val # # # La valeur de pi est passee par le common/const/ou /const2/ . # Sinon, il faut la calculer avant d'appeler ces fonctions . # # ---------------------------------------------------------------- # Fonctions a changer eventuellement, selon x(x) et y(y) choisis . # ----------------------------------------------------------------- # # ..... ici, on a l'application particuliere suivante ........ # # ************************************** # ** x = 2. * pi/iim * X ** # ** y = pi/jjm * Y ** # ************************************** # # .................................................................. # .................................................................. # # # #----------------------------------------------------------------------- def fxysinus (ddy,ddx): """c c Calcul des longitudes et des latitudes pour une fonction f(x,y) c avec y = Asin( j ) . c c Auteur : P. Le Van c c from: dyn3d/fxysinus.F """ fname = 'fxysinus' # ...... calcul des latitudes et de y' ..... # rlatu = np.zeros((ddy+1), dtype=np.float) yprimu = np.zeros((ddy+1), dtype=np.float) for j in range(ddy+1): rlatu[j] = fy( np.float(j+1), ddy) yprimu[j] = fyprim( np.float(j+1), ddy) rlatv = np.zeros((ddy), dtype=np.float) rlatu1 = np.zeros((ddy), dtype=np.float) rlatu2 = np.zeros((ddy), dtype=np.float) yprimv = np.zeros((ddy), dtype=np.float) yprimu1 = np.zeros((ddy), dtype=np.float) yprimu2 = np.zeros((ddy), dtype=np.float) for j in range(ddy): rlatv[j] = fy( np.float(j) + 0.5, ddy) rlatu1[j] = fy( np.float(j) + 0.25, ddy) rlatu2[j] = fy( np.float(j) + 0.75, ddy) yprimv[j] = fyprim( np.float(j) + 0.5, ddy) yprimu1[j] = fyprim( np.float(j) + 0.25, ddy) yprimu2[j] = fyprim( np.float(j) + 0.75, ddy) # # ..... calcul des longitudes et de x' ..... # rlonv = np.zeros((ddx+1), dtype=np.float) rlonu = np.zeros((ddx+1), dtype=np.float) rlonm025 = np.zeros((ddx+1), dtype=np.float) rlonp025 = np.zeros((ddx+1), dtype=np.float) xprimv = np.zeros((ddx+1), dtype=np.float) xprimu = np.zeros((ddx+1), dtype=np.float) xprimm025 = np.zeros((ddx+1), dtype=np.float) xprimp025 = np.zeros((ddx+1), dtype=np.float) for i in range(ddx + 1): rlonv[i] = fx( np.float(i), ddx ) rlonu[i] = fx( np.float(i)+ 0.5, ddx ) rlonm025[i] = fx( np.float(i)- 0.25, ddx ) rlonp025[i] = fx( np.float(i)+ 0.25, ddx ) xprimv[i] = fxprim(np.float(i), ddx ) xprimu[i] = fxprim(np.float(i)+ 0.5, ddx ) xprimm025[i] = fxprim(np.float(i)- 0.25, ddx ) xprimp025[i] = fxprim(np.float(i)+ 0.25, ddx ) return rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, \ xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 def fxhyp (dx, dy, xzoomdeg, grossism): """ c Auteur : P. Le Van c Calcule les longitudes et derivees dans la grille du GCM pour une c fonction f(x) a tangente hyperbolique . c c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.) c dzoom etant la distance totale de la zone du zoom c tau la raideur de la transition de l'interieur a l'exterieur du zoom c c On doit avoir grossism x dzoom < pi ( radians ) , en longitude. c ******************************************************************** """ fname = 'fxhyp' nmax = 30000 nmax2 = 2*nmax scal180 = True # scal180 = .TRUE. si on veut avoir le premier point scalaire pour # une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres. # sinon scal180 = .FALSE. # ...... arguments d'entree ....... # # REAL xzoomdeg,dzooma,tau,grossism # ...... arguments de sortie ...... rlonm025 = np.zeros((dx+1), dtype=np.float) xprimm025 = np.zeros((dx+1), dtype=np.float) rlonv = np.zeros((dx+1), dtype=np.float) xprimv = np.zeros((dx+1), dtype=np.float) rlonu = np.zeros((dx+1), dtype=np.float) xprimu = np.zeros((dx+1), dtype=np.float) rlonp025 = np.zeros((dx+1), dtype=np.float) xprimp025 = np.zeros((dx+1), dtype=np.float) # .... variables locales .... # xlon = np.zeros((dx+1), dtype=np.float32) xprimm = np.zeros((dx+1), dtype=np.float32) xtild = np.zeros((nmax2), dtype=np.float32) fhyp = np.zeros((nmax2), dtype=np.float32) Xprimt = np.zeros((nmax2), dtype=np.float32) Xf = np.zeros((nmax2), dtype=np.float32) xxpr = np.zeros((nmax2), dtype=np.float32) xvrai = np.zeros((dx+1), dtype=np.float32) xxprim = np.zeros((dx+1), dtype=np.float32) pi = np.pi depi = 2. * np.pi epsilon = 1.e-3 xzoom = xzoomdeg * pi/180. if dx == 1: rlonm025[1] = -pi/2. rlonv[1] = 0. rlonu[1] = pi rlonp025[1] = pi/2. rlonm025[2] = rlonm025[1] + depi rlonv[2] = rlonv[1] + depi rlonu[2] = rlonu[1] + depi rlonp025[2] = rlonp025[1] + depi xprimm025 = 1. xprimv = 1. xprimu = 1. xprimp025 = 1. champmin =depi champmax = depi return dzooma, tau, rlonm025, xprimm025, rlonv, xprimv, rlonu, xprimu, \ rlonp025, xprimp025, champmin, champmax decalx = .75 if grossism == 1. and scal180: decalx = 1. print 'FXHYP scal180,decalx', scal180,decalx if dzooma > 1.: dzoom = dzooma * depi elif dzooma < 25.: print erormsg print ' ' + fname +": Le param. dzoomx pour 'fxhyp' est trop petit dzooma:",\ dzooma,'!L augmenter et relancer !' quit(-1) else: dzoom = dzooma * pi/180. print ' xzoom( rad.),grossism,tau,dzoom (radians)' print xzoom,grossism,tau,dzoom xtild = - pi + np.range(namx2) * depi /nmax2 for i in range(nmax,nmax2): fa = tau* ( dzoom/2. - xtild[i] ) fb = xtild[i] * ( pi - xtild[i] ) if 200.* fb < - fa: fhyp[i] = -1. elif 200. * fb < fa: fhyp[i] = 1. else: if np.abs(fa) < 1.e-13 and np.abs(fb) < 1.e-13: if 200.*fb + fa < 1.e-10: fhyp[i] = -1. elif 200.*fb - fa < 1.e-10: fhyp[i] = 1. else: fhyp[i] = np.tanh(fa/fb) if xtild[i] == 0.: fhyp[i] = 1. if xtild[i] == pi: fhyp[i] = -1. ## .... Calcul de beta .... ffdx = 0. for i in range(nmax+1,nmax2): xmoy = 0.5 * ( xtild[i-1] + xtild[i] ) fa = tau*( dzoom/2. - xmoy ) fb = xmoy*( pi - xmoy ) if 200.* fb < -fa: fxm = -1. elif 200. * fb < fa: fxm = 1. else: if np.abs(fa) < 1.e-13 and np.abs(fb) < 1.e-13: if 200.*fb + fa < 1.e-10: fxm = -1. elif 200.*fb - fa < 1.e-10: fxm = 1. else: fxm = np.tanh( fa/fb ) if xmoy == 0.: fxm = 1. if xmoy == pi: fxm = -1. ffdx = ffdx + fxm * ( xtild[i] - xtild[i-1] ) beta = ( grossism * ffdx - pi ) / ( ffdx - pi ) if 2.*beta - grossism < 0.: print warnmsg print ' ' + fname + ': ** Attention ! La valeur beta calculee dans la ' + \ 'routine fxhyp est mauvaise ! ' print ' Modifier les valeurs de grossismx ,tau ou dzoomx et relancer ! ***' quit(-1) # # ..... calcul de Xprimt ..... # for i in range(nmax, nmax2): Xprimt[i] = beta + ( grossism - beta ) * fhyp[i] for i in range(nmax+1, nmax2): Xprimt[nmax2-i] = Xprimt[i] # ..... Calcul de Xf ........ Xf[0] = - pi for i in range(nmax+1, nmax2): xmoy = 0.5 * ( xtild[i-1] + xtild[i] ) fa = tau* ( dzoom/2. - xmoy ) fb = xmoy * ( pi - xmoy ) if 200.* fb < -fa: fxm = -1. elif 200. * fb < fa: fxm = 1. else: fxm = np.tanh( fa/fb ) if xmoy == 0.: fxm = 1. if xmoy == pi: fxm = -1. xxpr[i] = beta + ( grossism - beta ) * fxm for i in range(nmax+1, nmax2): xxpr[nmax2-i+1] = xxpr[i] for i in range(nmax2): Xf[i] = Xf[i-1] + xxpr[i] * ( xtild[i] - xtild[i-1] ) # ***************************************************************** # # ..... xuv = 0. si calcul aux pts scalaires ........ # ..... xuv = 0.5 si calcul aux pts U ........ for ik in range(4): if ik == 0: xuv = -0.25 elif ik == 1: xuv = 0. elif ik == 2: xuv = 0.50 elif ik == 4: xuv = 0.25 xo1 = 0. ii1=1 ii2=dx if ik == 0 and grossism == 1.: ii1 = 2 ii2 = dx+1 for i in range(ii1, ii2): xlon2 = - pi + (np.float(i) + xuv - decalx) * depi / np.float32(dx) Xfi = xlon2 ended = False for it in range(nmax2,0,-1): if Xfi > Xf[it]: ended = True break if not ended: it = 0 # ...... Calcul de Xf(xi) ...... # xi = xtild[it] if it == nmax2: it = nmax2 -1 Xf[it+1] = pi # ..................................................................... # # Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un # polynome de degre 3 qui passe par les points (Xf(it),xtild(it) ) # et (Xf(it+1),xtild(it+1) ) a0, a1, a2, a3 = coefpoly( Xf[it], Xf[it+1], Xprimt[it], Xprimt[it+1], \ xtild[it], xtild[it+1]) Xf1 = Xf[it] Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi for iteri in range(300): xi = xi - ( Xf1 - Xfi )/ Xprimin ended = False if np.abs(xi-xo1) < epsilon: ended = True break xo1 = xi xi2 = xi * xi Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2 if not ended: print errmsg print ' ' + fname + ': Pas de solution ***** ',i,xlon2,iteri quit(-1) xxprim[i] = depi/ ( np.float32(dx) * Xprimin ) xvrai[i] = xi + xzoom if ik == 0 and grossism == 1: xvrai[0] = xvrai[dx+1] - depi xxprim[0] = xxprim[dx+1] xlon[0:dx+1] = xvrai[0:dx+1] xprimm[0:dx+1] = xxprim[0:dx+1] for i in range(dx-1): if xvrai[i+1] < xvrai[i]: print errormsg print ' ' + fname + ': PBS. avec rlonu(',i+1,') plus petit que rlonu(', \ i,')' quit(-1) # # ... Reorganisation des longitudes pour les avoir entre - pi et pi .. # ........................................................................ champmin = 1.e12 champmax = -1.e12 for i in range(dx): champmin = np.min( [champmin,xvrai[i]] ) champmax = np.max( [champmax,xvrai[i]] ) # if champmin >= -pi-0.10 and champmax <= pi+0.10: # GO TO 1600 # # HERE -- here # # ELSE # WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi', # , ' et pi ' #c # IF( xzoom.LE.0.) THEN # IF( ik.EQ. 1 ) THEN # DO i = 1, iim # IF( xvrai(i).GE. - pi ) GO TO 80 # ENDDO # WRITE(6,*) ' PBS. 1 ! Xvrai plus petit que - pi ! ' # STOP 8 # 80 CONTINUE # is2 = i # ENDIF # # IF( is2.NE. 1 ) THEN # DO ii = is2 , iim # xlon (ii-is2+1) = xvrai(ii) # xprimm(ii-is2+1) = xxprim(ii) # ENDDO # DO ii = 1 , is2 -1 # xlon (ii+iim-is2+1) = xvrai(ii) + depi # xprimm(ii+iim-is2+1) = xxprim(ii) # ENDDO # ENDIF # ELSE # IF( ik.EQ.1 ) THEN # DO i = iim,1,-1 # IF( xvrai(i).LE. pi ) GO TO 90 # ENDDO # WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! ' # STOP 9 # 90 CONTINUE # is2 = i # ENDIF # idif = iim -is2 # DO ii = 1, is2 # xlon (ii+idif) = xvrai(ii) # xprimm(ii+idif) = xxprim(ii) # ENDDO # DO ii = 1, idif # xlon (ii) = xvrai (ii+is2) - depi # xprimm(ii) = xxprim(ii+is2) # ENDDO # ENDIF # ENDIF #c #c ......... Fin de la reorganisation ............................ # 1600 CONTINUE # xlon ( iip1) = xlon(1) + depi # xprimm( iip1 ) = xprimm (1 ) # DO i = 1, iim+1 # xvrai(i) = xlon(i)*180./pi # ENDDO # IF( ik.EQ.1 ) THEN #c WRITE(6,*) ' XLON aux pts. V-0.25 apres ( en deg. ) ' #c WRITE(6,18) #c WRITE(6,68) xvrai #c WRITE(6,*) ' XPRIM k ',ik #c WRITE(6,566) xprimm # DO i = 1,iim +1 # rlonm025(i) = xlon( i ) # xprimm025(i) = xprimm(i) # ENDDO # ELSE IF( ik.EQ.2 ) THEN #c WRITE(6,18) #c WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) ' #c WRITE(6,68) xvrai #c WRITE(6,*) ' XPRIM k ',ik #c WRITE(6,566) xprimm # DO i = 1,iim + 1 # rlonv(i) = xlon( i ) # xprimv(i) = xprimm(i) # ENDDO # ELSE IF( ik.EQ.3) THEN #c WRITE(6,18) #c WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) ' #c WRITE(6,68) xvrai #c WRITE(6,*) ' XPRIM ik ',ik #c WRITE(6,566) xprimm # DO i = 1,iim + 1 # rlonu(i) = xlon( i ) # xprimu(i) = xprimm(i) # ENDDO # ELSE IF( ik.EQ.4 ) THEN #c WRITE(6,18) #c WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) ' #c WRITE(6,68) xvrai #c WRITE(6,*) ' XPRIM ik ',ik #c WRITE(6,566) xprimm # DO i = 1,iim + 1 # rlonp025(i) = xlon( i ) # xprimp025(i) = xprimm(i) # ENDDO # ENDIF #5000 CONTINUE #c # WRITE(6,18) #c #c ........... fin de la boucle do 5000 ............ # DO i = 1, iim # xlon(i) = rlonv(i+1) - rlonv(i) # ENDDO # champmin = 1.e12 # champmax = -1.e12 # DO i = 1, iim # champmin = MIN( champmin, xlon(i) ) # champmax = MAX( champmax, xlon(i) ) # ENDDO # champmin = champmin * 180./pi # champmax = champmax * 180./pi #18 FORMAT(/) #24 FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3) #68 FORMAT(1x,7f9.2) #566 FORMAT(1x,7f9.4) return dzooma, tau, rlonm025, xprimm025, rlonv, xprimv, rlonu, xprimu, \ rlonp025, xprimp025, champmin, champmax def fxyhyper (ddy, ddx, yzoom, grossy, dzoomy, tauy, xzoom, grossx, dzoomx, taux, \ rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025): """c c Auteur : P. Le Van . c c d'apres formulations de R. Sadourny . c c c Ce spg calcule les latitudes( routine fyhyp ) et longitudes( fxhyp ) c par des fonctions a tangente hyperbolique . c c Il y a 3 parametres ,en plus des coordonnees du centre du zoom (xzoom c et yzoom ) : c c a) le grossissement du zoom : grossy ( en y ) et grossx ( en x ) c b) l' extension du zoom : dzoomy ( en y ) et dzoomx ( en x ) c c) la raideur de la transition du zoom : taux et tauy c c N.B : Il vaut mieux avoir : grossx * dzoomx < pi ( radians ) c ****** c et grossy * dzoomy < pi/2 ( radians ) c """ fname = 'fxyhyper' # CALL fyhyp ( yzoom, grossy, dzoomy,tauy , # , rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 , # , dymin,dymax ) # CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv, # , xprimv,rlonu,xprimu,rlonp025,xprimp025 , dxmin,dxmax ) for i in range(ddx+1): if rlonp025[i] < rlonv[i]: print errormsg print ' ' + fname + ' Attention ! rlonp025 < rlonv',i quit(-1) if rlonv[i] < rlonm025[i]: print errormsg print ' ' + fname + ' Attention ! rlonm025 > rlonv',i quit(-1) if rlonp025[i] > rlonu[i]: print errormsg print ' ' + fname + ' Attention ! rlonp025 > rlonu',i quit(-1) print ' *** TEST DE COHERENCE OK POUR FX **** ' for j in range(ddy): if rlatu1[j] <= rlatu2[j]: print errormsg print ' ' + fname + 'Attention ! rlatu1 < rlatu2',rlatu1[j],rlatu2[j],j quit(-1) if rlatu2[j] <= rlatu[j+1]: print errormsg print ' ' + fname + 'Attention ! rlatu2 < rlatup1',rlatu2[j],rlatu[j+1],j quit(-1) if rlatu[j] <= rlatu1[j]: print errormsg print ' ' + fname + 'Attention ! rlatu < rlatu1',rlatu[j],rlatu1[j],j quit(-1) if rlatv[j] <= rlatu2[j]: print errormsg print ' ' + fname + 'Attention ! rlatv < rlatu2 ',rlatv[j],rlatu2[j],j quit(-1) if rlatv[j] >= rlatu1[j]: print errormsg print ' ' + fname + 'Attention ! rlatv > rlatu1 ',rlatv[j],rlatu1[j],j quit(-1) if rlatv[j] >= rlatu[j]: print errormsg print ' ' + fname + 'Attention ! rlatv > rlatu ',rlatv[j],rlatu[j],j quit(-1) print ' *** TEST DE COHERENCE OK POUR FY **** ' print ' Latitudes ' print ' *********** ' print 'Au centre du zoom, la longueur de la maille est d environ', dymin , \ 'degres alors que la maille en dehors de la zone du zoom est d environ', \ dymax, 'degres' print ' Si cette derniere est trop lache , modifiez les parametres grossism , ' +\ 'tau , dzoom pour Y et repasser ! ' print ' Longitudes ' print ' ************ ' print 'Au centre du zoom, la longueur de la maille est d environ', dxmin , \ 'degres alors que la maille en dehors de la zone du zoom est d environ', \ dxmax, 'degres' print ' Si cette derniere est trop lache , modifiez les parametres grossism , ' +\ 'tau , dzoom pour Y et repasser ! ' return def inigeom(dy,dx): """c c Auteur : P. Le Van c c ............ Version du 01/04/2001 ........................ c c Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4 aux memes en- c endroits que les aires aireij1,..aireij4 . c Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol. c c """ fname = 'inigeom' cvu = np.zeros((dy+1,dx+1), dtype=np.float) cuv = np.zeros((dy, dx+1), dtype=np.float) cuij1 = np.zeros((dy+1,dx+1), dtype=np.float) cuij2 = np.zeros((dy+1,dx+1), dtype=np.float) cuij3 = np.zeros((dy+1,dx+1), dtype=np.float) cuij4 = np.zeros((dy+1,dx+1), dtype=np.float) cvij1 = np.zeros((dy+1,dx+1), dtype=np.float) cvij2 = np.zeros((dy+1,dx+1), dtype=np.float) cvij3 = np.zeros((dy+1,dx+1), dtype=np.float) cvij4 = np.zeros((dy+1,dx+1), dtype=np.float) aireij1 = np.zeros((dy+1,dx+1), dtype=np.float) aireij2 = np.zeros((dy+1,dx+1), dtype=np.float) aireij3 = np.zeros((dy+1,dx+1), dtype=np.float) aireij4 = np.zeros((dy+1,dx+1), dtype=np.float) aire = np.zeros((dy+1,dx+1), dtype=np.float) aireu = np.zeros((dy+1,dx+1), dtype=np.float) airev = np.zeros((dy+1,dx+1), dtype=np.float) unsaire = np.zeros((dy+1,dx+1), dtype=np.float) unsair_gam1 = np.zeros((dy+1,dx+1), dtype=np.float) unsair_gam2 = np.zeros((dy+1,dx+1), dtype=np.float) airesurg = np.zeros((dy+1,dx+1), dtype=np.float) unsairez = np.zeros((dy+1,dx+1), dtype=np.float) unsairz_gam = np.zeros((dy+1,dx+1), dtype=np.float) fext = np.zeros((dy+1,dx+1), dtype=np.float) alpha1 = np.zeros((dy+1,dx+1), dtype=np.float) alpha2 = np.zeros((dy+1,dx+1), dtype=np.float) alpha3 = np.zeros((dy+1,dx+1), dtype=np.float) alpha4 = np.zeros((dy+1,dx+1), dtype=np.float) alpha1p2 = np.zeros((dy+1,dx+1), dtype=np.float) alpha1p4 = np.zeros((dy+1,dx+1), dtype=np.float) alpha2p3 = np.zeros((dy+1,dx+1), dtype=np.float) alpha3p4 = np.zeros((dy+1,dx+1), dtype=np.float) rlonvv = np.zeros((dx+1), dtype=np.float) rlatuu = np.zeros((dy+1), dtype=np.float) rlatu1 = np.zeros((dy), dtype=np.float) yprimu1 = np.zeros((dy), dtype=np.float) rlatu2 = np.zeros((dy), dtype=np.float) yprimu2 = np.zeros((dy), dtype=np.float) yprimv = np.zeros((dy), dtype=np.float) yprimu = np.zeros((dy+1), dtype=np.float) rlonm025 = np.zeros((dx+1), dtype=np.float) xprimm025 = np.zeros((dx+1), dtype=np.float) rlonp025 = np.zeros((dx+1), dtype=np.float) xprimp025 = np.zeros((dx+1), dtype=np.float) cu = np.zeros((dy+1,dx+1), dtype=np.float) cv = np.zeros((dy+1,dx+1), dtype=np.float) unscu2 = np.zeros((dy+1,dx+1), dtype=np.float) unscv2 = np.zeros((dy+1,dx+1), dtype=np.float) cuvsurcv = np.zeros((dy+1,dx+1), dtype=np.float) cuvsurcuv = np.zeros((dy+1,dx+1), dtype=np.float) cvsurcv = np.zeros((dy+1,dx+1), dtype=np.float) cvsurcuv = np.zeros((dy+1,dx+1), dtype=np.float) cuvscvgam1 = np.zeros((dy+1,dx+1), dtype=np.float) cuvscvgam2 = np.zeros((dy+1,dx+1), dtype=np.float) cvscuvgam = np.zeros((dy+1,dx+1), dtype=np.float) cvusurcu = np.zeros((dy+1,dx+1), dtype=np.float) cusurcvu = np.zeros((dy+1,dx+1), dtype=np.float) cvuscugam1 = np.zeros((dy+1,dx+1), dtype=np.float) cvuscugam2 = np.zeros((dy+1,dx+1), dtype=np.float) cuscvugam = np.zeros((dy+1,dx+1), dtype=np.float) airvscu2 = np.zeros((dy+1,dx+1), dtype=np.float) aivscu2gam = np.zeros((dy+1,dx+1), dtype=np.float) airuscv2 = np.zeros((dy+1,dx+1), dtype=np.float) aiuscv2gam = np.zeros((dy+1,dx+1), dtype=np.float) constang = np.zeros((dy+1,dx+1), dtype=np.float) # # # ------------------------------------------------------------------ # - - # - calcul des coeff. ( cu, cv , 1./cu**2, 1./cv**2 ) - # - - # ------------------------------------------------------------------ # # les coef. ( cu, cv ) permettent de passer des vitesses naturelles # aux vitesses covariantes et contravariantes , ou vice-versa ... # # # on a : u (covariant) = cu * u (naturel) , u(contrav)= u(nat)/cu # v (covariant) = cv * v (naturel) , v(contrav)= v(nat)/cv # # on en tire : u(covariant) = cu * cu * u(contravariant) # v(covariant) = cv * cv * v(contravariant) # # # on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2 # = = # et - jm/2 < Y < jm/2 # = = # # ................................................... # ................................................... # . x est la longitude du point en radians . # . y est la latitude du point en radians . # . . # . on a : cu(i,j) = rad * COS(y) * dx/dX . # . cv( j ) = rad * dy/dY . # . aire(i,j) = cu(i,j) * cv(j) . # . . # . y, dx/dX, dy/dY calcules aux points concernes . # . . # ................................................... # ................................................... # # # # , # cv , bien que dependant de j uniquement,sera ici indice aussi en i # pour un adressage plus facile en ij . # # # # ************** aux points u et v , ***************** # xprimu et xprimv sont respectivement les valeurs de dx/dX # yprimu et yprimv . . . . . . . . . . . dy/dY # rlatu et rlatv . . . . . . . . . . .la latitude # cvu et cv . . . . . . . . . . . cv # # ************** aux points u, v, scalaires, et z **************** # cu, cuv, cuscal, cuz sont respectiv. les valeurs de cu # # # # Exemple de distribution de variables sur la grille dans le # domaine de travail ( X,Y ) . # ................................................................ # DX=DY= 1 # # # + represente un point scalaire ( p.exp la pression ) # > represente la composante zonale du vent # V represente la composante meridienne du vent # o represente la vorticite # # ---- , car aux poles , les comp.zonales covariantes sont nulles # # # # i -> # # 1 2 3 4 5 6 7 8 # j # v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + -- # # V o V o V o V o V o V o V o V o # # 2 + > + > + > + > + > + > + > + > # # V o V o V o V o V o V o V o V o # # 3 + > + > + > + > + > + > + > + > # # V o V o V o V o V o V o V o V o # # 4 + > + > + > + > + > + > + > + > # # V o V o V o V o V o V o V o V o # # 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + -- # # # Ci-dessus, on voit que le nombre de pts.en longitude est egal # a IM = 8 # De meme , le nombre d'intervalles entre les 2 poles est egal # a JM = 4 # # Les points scalaires ( + ) correspondent donc a des valeurs # entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) . # # Les vents U ( > ) correspondent a des valeurs semi- # entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1) # # Les vents V ( V ) correspondent a des valeurs entieres # de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5) # # # if nitergdiv != 2: gamdi_gdiv = coefdis/ ( np.float(nitergdiv) -2. ) else: gamdi_gdiv = 0. if nitergrot != 2: gamdi_grot = coefdis/ ( np.float(nitergrot) -2. ) else: gamdi_grot = 0. if niterh != 2: gamdi_h = coefdis/ ( np.float(niterh) -2. ) else: gamdi_h = 0. print 'gamdi_gd:',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,nitergdiv,nitergrot,niterh # ---------------------------------------------------------------- # if not fxyhypb: if ysinus: print ' *** Inigeom , Y = Sinus ( Latitude ) *** ' # .... utilisation de f(x,y ) avec y = sinus de la latitude ..... rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, \ xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 = \ fxysinus(dx, dy) else: print '*** Inigeom , Y = Latitude , der. sinusoid . ***' # .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ... pxo = clon *np.pi /180. pyo = 2.* clat* np.pi /180. # # .... determination de transx ( pour le zoom ) par Newton-Raphson ... # itmax = 10 eps = .1e-7 xo1 = 0. for iter in range(itmax): x1 = xo1 f = x1+ alphax*np.sin(x1-pxo) df = 1.+ alphax*np.cos(x1-pxo) x1 = x1 - f/df xdm = np.abs( x1- xo1 ) if xdm > eps: xo1 = x1 transx = xo1 itmay = 10 eps = .1e-7 yo1 = 0. for iter in range(itmay): y1 = yo1 f = y1 + alphay*np.sin(y1-pyo) df = 1. + alphay*np.cos(y1-pyo) y1 = y1 -f/df ydm = np.abs(y1-yo1) if ydm > eps: yo1 = y1 transy = yo1 rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, \ xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 = \ fxy(dx,dy) else: # # .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol. # ..................................................................... print '*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' #!!!!! Lluis #!! HERE, how to do all this without zoom!? #! # CALL fxyhyper( clat, grossismy, dzoomy, tauy , # , clon, grossismx, dzoomx, taux , # , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 , # , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 ) # ------------------------------------------------------------------- rlatu[0] = np.arcsin(1.) rlatu[dy] = -rlatu[0] # .... calcul aux poles .... yprimu[0] = 0. yprimu[dy] = 0. un4rad2 = 0.25 * rad * rad # # -------------------------------------------------------------------- # -------------------------------------------------------------------- # - - # - calcul des aires ( aire,aireu,airev, 1./aire, 1./airez ) - # - et de fext , force de coriolis extensive . - # - - # -------------------------------------------------------------------- # -------------------------------------------------------------------- # # # # A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont # affectees 4 aires entourant P , calculees respectivement aux points # ( i + 1/4, j - 1/4 ) : aireij1 (i,j) # ( i + 1/4, j + 1/4 ) : aireij2 (i,j) # ( i - 1/4, j + 1/4 ) : aireij3 (i,j) # ( i - 1/4, j - 1/4 ) : aireij4 (i,j) # # , # Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y). # Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme # des 4 aires aireij1,aireij2,aireij3,aireij4 qui sont affectees au # point (i,j) . # On definit en outre les coefficients alpha comme etant egaux a # (aireij / aire), c.a.d par exp. alpha1(i,j)=aireij1(i,j)/aire(i,j) # # De meme, toute aire centree en 1 point U est egale a la somme des # 4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U . # Idem pour airev, airez . # # On a ,pour chaque maille : dX = dY = 1 # # # . V # # aireij4 . . aireij1 # # U . . P . U # # aireij3 . . aireij2 # # . V # # # # # # .................................................................... # # Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4 # qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen # taires cuij et les 4 elongat. cvij qui sont calculees aux memes # endroits que les aireij . # # .................................................................... # # ....... do 35 : boucle sur les jjm + 1 latitudes ..... # # for j in range(dy+1): if j == 1: yprm = yprimu1[j] rlatm = rlatu1[j] coslatm = np.cos( rlatm ) radclatm = 0.5* rad * coslatm for i in range(dx): xprp = xprimp025[i] xprm = xprimm025[i] aireij2[0,i] = un4rad2 * coslatm * xprp * yprm aireij3[0,i] = un4rad2 * coslatm * xprm * yprm cuij2[0,i] = radclatm * xprp cuij3[0,i] = radclatm * xprm cvij2[0,i] = 0.5* rad * yprm cvij3[0,i] = cvij2[0,i] for i in range(dx): aireij1[0,i] = 0. aireij4[0,i] = 0. cuij1[0,i] = 0. cuij4[0,i] = 0. cvij1[0,i] = 0. cvij4[0,i] = 0. elif j == dy: yprp = yprimu2[j-1] rlatp = rlatu2[j-1] coslatp = np.cos( rlatp ) radclatp = 0.5* rad * coslatp for i in range(dx): xprp = xprimp025[i] xprm = xprimm025[i] aireij1[dy,i] = un4rad2 * coslatp * xprp * yprp aireij4[dy,i] = un4rad2 * coslatp * xprm * yprp cuij1[dy,i] = radclatp * xprp cuij4[dy,i] = radclatp * xprm cvij1[dy,i] = 0.5 * rad* yprp cvij4[dy,i] = cvij1[dy,i] for i in range(dx): aireij2[dy,i] = 0. aireij3[dy,i] = 0. cvij2[dy,i] = 0. cvij3[dy,i] = 0. cuij2[dy,i] = 0. cuij3[dy,i] = 0. else: rlatp = rlatu2[j-1] yprp = yprimu2[j-1] rlatm = rlatu1[j] yprm = yprimu1[j] coslatm = np.cos( rlatm ) coslatp = np.cos( rlatp ) radclatp = 0.5* rad * coslatp radclatm = 0.5* rad * coslatm for i in range(dx): xprp = xprimp025[i] xprm = xprimm025[i] ai14 = un4rad2 * coslatp * yprp ai23 = un4rad2 * coslatm * yprm aireij1[j,i] = ai14 * xprp aireij2[j,i] = ai23 * xprp aireij3[j,i] = ai23 * xprm aireij4[j,i] = ai14 * xprm cuij1[j,i] = radclatp * xprp cuij2[j,i] = radclatm * xprp cuij3[j,i] = radclatm * xprm cuij4[j,i] = radclatp * xprm cvij1[j,i] = 0.5* rad * yprp cvij2[j,i] = 0.5* rad * yprm cvij3[j,i] = cvij2[j,i] cvij4[j,i] = cvij1[j,i] # # ........ periodicite ............ # cvij1[j,dx] = cvij1[j,0] cvij2[j,dx] = cvij2[j,0] cvij3[j,dx] = cvij3[j,0] cvij4[j,dx] = cvij4[j,0] cuij1[j,dx] = cuij1[j,0] cuij2[j,dx] = cuij2[j,0] cuij3[j,dx] = cuij3[j,0] cuij4[j,dx] = cuij4[j,0] aireij1[j,dx] = aireij1[j,0] aireij2[j,dx] = aireij2[j,0] aireij3[j,dx] = aireij3[j,0] aireij4[j,dx] = aireij4[j,0] # # .............................................................. # for j in range(dy+1): for i in range(dx): aire[j,i] = aireij1[j,i] + aireij2[j,i] + aireij3[j,i] + aireij4[j,i] alpha1[j,i] = aireij1[j,i] / aire[j,i] alpha2[j,i] = aireij2[j,i] / aire[j,i] alpha3[j,i] = aireij3[j,i] / aire[j,i] alpha4[j,i] = aireij4[j,i] / aire[j,i] alpha1p2[j,i] = alpha1 [j,i] + alpha2 [j,i] alpha1p4[j,i] = alpha1 [j,i] + alpha4 [j,i] alpha2p3[j,i] = alpha2 [j,i] + alpha3 [j,i] alpha3p4[j,i] = alpha3 [j,i] + alpha4 [j,i] aire[j,dx] = aire[j,0] alpha1[j,dx] = alpha1[j,0] alpha2[j,dx] = alpha2[j,0] alpha3[j,dx] = alpha3[j,0] alpha4[j,dx] = alpha4[j,0] alpha1p2[j,dx] = alpha1p2[j,0] alpha1p4[j,dx] = alpha1p4[j,0] alpha2p3[j,dx] = alpha2p3[j,0] alpha3p4[j,dx] = alpha3p4[j,0] for j in range(dy+1): for i in range(dx): aireu[j,i] = aireij1[j,i] + aireij2[j,i] + aireij4[j,i+1] + aireij3[j,i+1] unsaire[j,i] = 1./ aire[j,i] unsair_gam1[j,i] = unsaire[j,i]** ( - gamdi_gdiv ) unsair_gam2[j,i] = unsaire[j,i]** ( - gamdi_h ) airesurg[j,i] = aire[j,i]/ g aireu[j,dx] = aireu[j,0] unsaire[j,dx] = unsaire[j,0] unsair_gam1[j,dx] = unsair_gam1[j,0] unsair_gam2[j,dx] = unsair_gam2[j,0] airesurg[j,dx] = airesurg[j,0] for j in range(dy): for i in range(dx): airev[j,i] = aireij2[j,i]+ aireij3[j,i]+ aireij1[j,i] + aireij4[j+1,i] for i in range(dx): airez = aireij2[j,i]+aireij1[j+1,i]+aireij3[j,i+1] + aireij4[j+1,i+1] unsairez[j,i] = 1./ airez unsairz_gam[j,i]= unsairez[j,i]** ( - gamdi_grot ) fext[j,i] = airez * np.sin(rlatv[j])* 2.* omeg airev[j,dx] = airev[j,0] unsairez[j,dx] = unsairez[j,0] fext[j,dx] = fext[j,0] unsairz_gam[j,dx] = unsairz_gam[j,0] # # ..... Calcul des elongations cu,cv, cvu ......... # for j in range(dy): for i in range(dx): cv[j,i] = 0.5 *( cvij2[j,i]+cvij3[j,i]+cvij1[j+1,i]+cvij4[j+1,i]) cvu[j,i]= 0.5 *( cvij1[j,i]+cvij4[j,i]+cvij2[j,i]+cvij3[j,i] ) cuv[j,i]= 0.5 *( cuij2[j,i]+cuij3[j,i]+cuij1[j+1,i]+cuij4[j+1,i]) unscv2[j,i] = 1./ ( cv[j,i]*cv[j,i] ) for i in range(dx): cuvsurcv [j,i] = airev[j,i] * unscv2[j,i] cvsurcuv [j,i] = 1./cuvsurcv[j,i] cuvscvgam1[j,i] = cuvsurcv [j,i] ** ( - gamdi_gdiv ) cuvscvgam2[j,i] = cuvsurcv [j,i] ** ( - gamdi_h ) cvscuvgam[j,i] = cvsurcuv [j,i] ** ( - gamdi_grot ) cv[j,dx] = cv[j,0] cvu[j,dx] = cvu[j,0] unscv2[j,dx] = unscv2[j,0] cuv[j,dx] = cuv[j,0] cuvsurcv[j,dx] = cuvsurcv[j,0] cvsurcuv[j,dx] = cvsurcuv[j,0] cuvscvgam1[j,dx] = cuvscvgam1[j,0] cuvscvgam2[j,dx] = cuvscvgam2[j,0] cvscuvgam[j,dx] = cvscuvgam[j,0] for j in range(1,dy): for i in range(dx): cu[j,i] = 0.5*(cuij1[j,i]+cuij4[j,i+1]+cuij2[j,i]+cuij3[j,i+1]) unscu2[j,i] = 1./ ( cu[j,i] * cu[j,i] ) cvusurcu[j,i] = aireu[j,i] * unscu2[j,i] cusurcvu[j,i] = 1./ cvusurcu[j,i] cvuscugam1[j,i] = cvusurcu[j,i] ** ( - gamdi_gdiv ) cvuscugam2[j,i] = cvusurcu[j,i] ** ( - gamdi_h ) cuscvugam[j,i] = cusurcvu[j,i] ** ( - gamdi_grot ) cu[j,dx] = cu[j,0] unscu2[j,dx] = unscu2[j,0] cvusurcu[j,dx] = cvusurcu[j,0] cusurcvu[j,dx] = cusurcvu[j,0] cvuscugam1[j,dx] = cvuscugam1[j,0] cvuscugam2[j,dx] = cvuscugam2[j,0] cuscvugam[j,dx] = cuscvugam[j,0] # # .... calcul aux poles .... # for i in range(dx+1): cu[0, i] = 0. unscu2[0, i] = 0. cvu[0, i] = 0. cu[dy,i] = 0. unscu2[dy,i] = 0. cvu[dy,i] = 0. # # .............................................................. # for j in range(dy): for i in range(dx): airvscu2[j,i] = airev[j,i]/ ( cuv[j,i] * cuv[j,i] ) aivscu2gam[j,i] = airvscu2[j,i]** ( - gamdi_grot ) airvscu2[j,dx] = airvscu2[j,0] aivscu2gam[j,dx] = aivscu2gam[j,0] for j in range(dy): for i in range(dx): airuscv2[j,i] = aireu[j,i]/ ( cvu[j,i] * cvu[j,i] ) aiuscv2gam[j,i] = airuscv2[j,i]** ( - gamdi_grot ) airuscv2[j,dx] = airuscv2[j,0] aiuscv2gam[j,dx] = aiuscv2gam[j,0] # # calcul des aires aux poles : # ----------------------------- # apoln = SSUM(dx,aire[0,0],1) apols = SSUM(dx,aire[dy,0],1) unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) ) unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) ) unsapolnga2 = 1./ ( apoln ** ( - gamdi_h ) ) unsapolsga2 = 1./ ( apols ** ( - gamdi_h ) ) # #----------------------------------------------------------------------- # gtitre='Coriolis version ancienne' # gfichier='fext1' # CALL writestd(fext,iip1*jjm) # # changement F. Hourdin calcul conservatif pour fext # constang contient le produit a * cos ( latitude ) * omega # for i in range(dx): constang[0,i] = 0. for j in range(dy-1): for i in range(dx): constang[j+1,i] = rad*omeg*cu[j+1,i]*np.cos(rlatu[j+1]) for i in range(dx): constang[dy,i] = 0. # # periodicite en longitude # for j in range(dy): fext[j,dx] = fext[j,0] for j in range(dy+1): constang[j,dx] = constang[j,0] # fin du changement # #----------------------------------------------------------------------- # print ' *** Coordonnees de la grille *** ' print ' LONGITUDES aux pts. V ( degres ) ' for i in range(dx+1): rlonvv[i] = rlonv[i]*180./np.pi print rlonvv print ' LATITUDES aux pts. V ( degres ) ' for i in range(dy): rlatuu[i] = rlatv[i]*180./np.pi print rlatuu[dy] for i in range(dx+1): rlonvv[i]=rlonu[i]*180./np.pi print ' LONGITUDES aux pts. U ( degres ) ' print rlonvv print ' LATITUDES aux pts. U ( degres ) ' for i in range(dy+1): rlatuu[i]=rlatu[i]*180./np.pi print rlatuu[0:dy+2] return aire, apoln, apols, airesurg, rlatu, rlatv, cu, cv def SSUM(n,sx,incx): """ Obsolete version of sum for non Fortan 90 code from dyn3d/cray.F """ ssumv = 0. ix = 0 if len(sx.shape) != 0: for i in range(n): ssumv=ssumv+sx[ix] ix=ix+incx else: ssumv=ssumv+sx return ssumv def SCOPY(ddx,ddy,ddz,sx,incx,sy,incy): """ Obsolete function to copy matrix values from dyn3d/cray.F """ fname = 'SCOPY' if len(sx.shape) == 2: for j in range(ddy-1): for i in range(ddx-1): sy[iy,ix] = sx[iy,ix] ix = ix+incx iy = iy+incy elif len(sx.shape) == 3: iy = 0 for j in range(ddy): ix = 0 for i in range(ddx): for l in range(ddz): sy[l,iy,ix] = sx[l,iy,ix] ix = ix+incx iy = iy+incy return sy def exner_hyb (dx, dy, dz, psv, pv, aire, apoln, apols): """c c Auteurs : P.Le Van , Fr. Hourdin . c .......... c c .... ngrid, ps,p sont des argum.d'entree au sous-prog ... c .... alpha,beta, pks,pk,pkf sont des argum.de sortie au sous-prog ... c c ************************************************************************ c Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des c couches . Pk(l) sera calcule aux milieux des couches l ,entre les c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . c ************************************************************************ c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont c la pression et la fonction d'Exner au sol . c c -------- z c A partir des relations ( 1 ) p*dz(pk) = kappa *pk*dz(p) et c ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1) c ( voir note de Fr.Hourdin ) , c c on determine successivement , du haut vers le bas des couches, les c coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), c puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, c pk(ij,l) donne par la relation (2), pour l = 2 a l = llm . c """ fname = 'exner_hyb' pksv = np.zeros((dy+1, dx+1), dtype=np.float) pkv = np.zeros((dz, dy+1, dx+1), dtype=np.float) pkfv = np.zeros((dz, dy+1, dx+1), dtype=np.float) ppn = np.zeros((dy+1, dx+1), dtype=np.float) pps = np.zeros((dy+1, dx+1), dtype=np.float) alphav = np.zeros((dz+1, dy+1, dx+1), dtype=np.float) betav = np.zeros((dz+1, dy+1, dx+1), dtype=np.float) if dz == 1: # Compute pks(:),pk(:),pkf(:) pks = (cpp/preff)*ps pk[0,:,:] = 0.5*pks # CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) is the same as the next line? pkf = pk # No filtering... not necessary on aquaplanet # CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) # our work is done, exit routine return pksv, pkv, pkfv #### General case: unpl2k = 1.+ 2.* kappa for j in range(dy+1): for i in range(dx+1): pksv[j,i] = cpp * ( psv[j,i]/preff ) ** kappa for j in range(dy+1): for i in range(dx+1): ppn[j,i] = aire[j,i] * pksv[j,i] pps[j,i] = aire[dy,i] * pksv[dy,i] xpn = SSUM(dx,ppn,1) /apoln xps = SSUM(dx,pps,1) /apols for j in range(dy): for i in range(dx): pksv[j,i] = xpn[i] pksv[dy-1,i] = xps[i] # # # .... Calcul des coeff. alpha et beta pour la couche l = llm .. # for j in range(dy+1): for i in range(dx+1): alphav[dz-1,j,i] = 0. betav[dz-1,j,i] = 1./ unpl2k # # ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ... # for l in range (dz-1,1,-1): for j in range(dy+1): for i in range(dx+1): dellta = pv[l,j,i]* unpl2k + pv[l+1,j,i]* ( betav[l+1,j,i]-unpl2k ) alphav[l,j,i] = -pv[l+1,j,i] / dellta * alphav[l+1,j,i] betav[l,j,i] = pv[l,j,i] / dellta # *********************************************************************** # ..... Calcul de pk pour la couche 1 , pres du sol .... # for j in range(dy+1): for i in range(dx+1): pkv[0,j,i] = ( pv[0,j,i]*pksv[j,i] - 0.5*alphav[1,j,i]*pv[1,j,i] ) \ *( pv[0,j,i]* (1.+kappa) + 0.5*( betav[1,j,i]-unpl2k )* pv[1,j,i] ) # # ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........ # for l in range(dz): for j in range(dy+1): for i in range(dx+1): pkv[l,j,i] = alphav[l,j,i] + betav[l,j,i] * pkv[l-1,j,i] # # pkfv = SCOPY ( dx+1, dy+1, dz, pkv, 1, pkfv, 1 ) # We do not filter for iniaqua # CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) return pksv, pkv, pkfv, alphav, betav def exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ): """c c Auteurs : F. Forget , Y. Wanherdrick c P.Le Van , Fr. Hourdin . c .......... c c .... ngrid, ps,p sont des argum.d'entree au sous-prog ... c .... beta, pks,pk,pkf sont des argum.de sortie au sous-prog ... c c ************************************************************************ c Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des c couches . Pk(l) sera calcule aux milieux des couches l ,entre les c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches . c ************************************************************************ c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont c la pression et la fonction d'Exner au sol . c c WARNING : CECI est une version speciale de exner_hyb originale c Utilise dans la version martienne pour pouvoir c tourner avec des coordonnees verticales complexe c => Il ne verifie PAS la condition la proportionalite en c energie totale/ interne / potentielle (F.Forget 2001) c ( voir note de Fr.Hourdin ) , c """ fname = 'exner_milieu' pks = np.zeros((dy, dx), dtype=np.float) pk = np.zeros((dz, dy, dx), dtype=np.float) pkf = np.zeros((dz, dy, dx), dtype=np.float) ppn = np.zeros((iim), dtype=np.float) ppn = np.zeros((iim), dtype=np.float) ip1jm = (dx+1)*dy firstcall = True modname = 'exner_milieu' # Sanity check if firstcall: # sanity checks for Shallow Water case (1 vertical layer) if llm == 1: if kappa != 1: print errormsg print ' ' + fname+ ': kappa!=1 , but running in Shallow Water mode!!' quit(-1) if cpp != r: print errormsg print ' ' + fname+ ': cpp!=r , but running in Shallow Water mode!!' quit(-1) firstcall = False #### Specific behaviour for Shallow Water (1 vertical layer) case: if llm == 1: # Compute pks(:),pk(:),pkf(:) for j,i in range(ngrid): pks[j,i] = (cpp/preff) * ps[j,i] pk[j,i,1] = .5*pks[j,i] pkf = SCOPY(dx,dy,dz, pk, 1, pkf, 1 ) # We do not filter for iniaqua # CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) # our work is done, exit routine return pksv, pkv, pkfv #### General case: # ------------- # Calcul de pks # ------------- for j,i in range(ngrid): pks[j,i] = cpp * ( ps[j,i]/preff ) ** kappa for j,i in range(iim): ppn[j,i] = aire[j,i] * pks[j,i] pps[j,i] = aire[j,i+ip1jm] * pks[j,i+ip1jm] xpn = SSUM(iim,ppn,1) /apoln xps = SSUM(iim,pps,1) /apols for j,i in range(iip1): pks[j,i] = xpn pks[j,i+ip1jm] = xps # # # .... Calcul de pk pour la couche l # -------------------------------------------- # dum1 = cpp * (2*preff)**(-kappa) for l in range(llm-1): for j,i in range(ngrid): pk[j,i,l] = dum1 * (p[j,i,l] + p[j,i,l+1])**kappa # .... Calcul de pk pour la couche l = llm .. # (on met la meme distance (en log pression) entre Pk(llm) # et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2) for j,i in range(ngrid): pk[j,i,llm] = pk[j,i,llm-1]**2 / pk[j,i,llm-2] # calcul de pkf # ------------- pkf = SCOPY( dx,dy,dz, pk, 1, pkf, 1 ) # We do not filter iniaqua # CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) # EST-CE UTILE ?? : calcul de beta # -------------------------------- for l in range(1, llm): for j,i in range(ngrid): beta[j,i,l] = pk[j,i,l] / pk[j,i,l-1] return pksv, pkv, pkfv def pression(dx, dy, dz, apv, bpv, psv): """c c Auteurs : P. Le Van , Fr.Hourdin . c ************************************************************************ c Calcule la pression p(l) aux differents niveaux l = 1 ( niveau du c sol) a l = llm +1 ,ces niveaux correspondant aux interfaces des (llm) c couches , avec p(ij,llm +1) = 0. et p(ij,1) = ps(ij) . c ************************************************************************ c """ fname = 'pression' press = np.zeros((dz+1, dy+1, dx+1), dtype=np.float) for l in range(dz+1): press[l,:,:] = apv[l] + bpv[l]*psv return press def sig_hybrid(sig,pa,preff): """ Function utilisee pour calculer des valeurs de sigma modifie pour conserver les coordonnees verticales decrites dans esasig.def/z2sig.def lors du passage en coordonnees hybrides F. Forget 2002 sig= sigma level pa= preff= reference pressure Connaissant sig (niveaux "sigma" ou on veut mettre les couches) L'objectif est de calculer newsig telle que (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig Cela ne se resoud pas analytiquement: => on resoud par iterration bourrine ---------------------------------------------- Information : where exp(1-1./x**2) become << x x exp(1-1./x**2) /x 1 1 0.68 0.5 0.5 1.E-1 0.391 1.E-2 0.333 1.E-3 0.295 1.E-4 0.269 1.E-5 0.248 1.E-6 => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25 """ fname = 'sig_hybrid' # maximum number of iterations maxiter = 9999 nsig = sig x1=0 x2=1 if sig >= 1.: nsig = sig elif sig*preff/pa >= 0.25: for j in range(maxiter): F = ((1. -pa/preff)*np.exp(1.-1./nsig**2)+(pa/preff)*nsig)/sig # print J,'nsig=', newsig, 'F=', F if F > 1: X2 = newsig nsig = (X1+nsig)*0.5 else: X1 = newsig nsig = (X2+nsig)*0.5 # Test : on arete lorsque on approxime sig a moins de 0.01 m pres # (en pseudo altitude) : if np.abs(10.*np.log(F)) < 1.e-5: break else: nsig= sig*preff/pa return nsig def presnivs_calc(vert_sampling, dz): """ From dyn3d/disvert.F calculation of vertical pressure levels vert_sampling= which kind of vertical sampling is desired: "param", "tropo", "strato" and "read" dz= numbef of vertical layers """ fname = 'presnivs_calc' llmp1 = dz + 1 pnivs = np.zeros((dz), dtype=np.float) s = np.zeros((dz), dtype=np.float) sig = np.zeros((dz+1), dtype=np.float) dsig = np.zeros((dz), dtype=np.float) dpres = np.zeros((dz), dtype=np.float) apv = np.zeros((dz+1), dtype=np.float) bpv = np.zeros((dz+1), dtype=np.float) # default scaleheight is 8km for earth scaleheight = 8. # vert_sampling = merge("strato", "tropo ", ok_strato) ! default value if vert_sampling == 'param': # On lit les options dans sigma.def: if not os.path.isfile('easig.def'): print errormsg print ' ' + fname + ": parameters file 'easig.def' does not exist!!" quit(-1) sfobj = open('sigma.def', 'r') scaleheight = np.float( ncvar.reduce_spaces(fobj.readline())) deltaz = np.float( ncvar.reduce_spaces(fobj.readline())) beta = np.float( ncvar.reduce_spaces(fobj.readline())) k0 = np.float( ncvar.reduce_spaces(fobj.readline())) k1 = np.float( ncvar.reduce_spaces(fobj.readline())) sfobj.close() alpha = deltaz/(dz*scaleheight) print ':scaleheight, alpha, k0, k1, beta', scaleheight, alpha, k0, k1, beta alpha=deltaz/np.tanh(1./k0)*2. zkm1=0. sig[0]=1. for l in range(dz): sig[l+1]=(np.cosh(l/k0))**(-alpha*k0/scaleheight) \ *exp(-alpha/scaleheight*np.tanh((llm-k1)/k0) \ *beta**(l-(llm-k1))/np.log(beta)) zk=-scaleheight*np.log(sig[l+1]) dzk1=alpha*np.tanh(l/k0) dzk2=alpha*np.tanh((llm-k1)/k0)*beta**(l-(llm-k1))/np.log(beta) print l, sig(l+1), zk, zk-zkm1, dzk1, dzk2 zkm1=zk sig[dz-1]=0. bpv[0:dz] = np.exp(1.-1./sig[0:dz]**2) bpv[llmp1-1] = 0. apv = pa * (sig - bp) elif vert_sampling == 'tropo': for l in range(dz): x = 2.*np.arcsin(1.)*(l-0.5)/(dz+1.) dsig[l] = 1.0+7.0*np.sin(x)**2 dsig = dsig/np.sum(dsig) sig[dz] = 0. for l in range(dz-1,0,-1): sig[l] = sig[l+1] + dsig[l] bpv[0]=1. bpv[1:dz] = np.exp(1.-1./sig[1:dz]**2) bpv[llmp1-1] = 0. apv[0] = 0. apv[1:dz+1] = pa*(sig[1:dz+1]-bpv[1:dz+1]) elif vert_sampling == 'strato': if dz == 39: dsigmin = 0.3 elif dz == 50: dsigmin = 1. else: print ' ATTENTION discretisation z a ajuster' dsigmin = 1. print 'Discretisation verticale DSIGMIN=',dsigmin for l in range(dz): x = 2.*np.arcsin(1.)*(l - 0.5)/(dz+1) dsig[l] =(dsigmin+7.*np.sin(x)**2) \ *(0.5*(1.-np.tanh(1.*(x-np.arcsin(1.))/np.arcsin(1.))))**2 dsig = dsig/np.sum(dsig) sig[dz] = 0. for l in range(dz-1,0,-1): sig[l] = sig[l+1] + dsig[l] bpv[0] = 1. bpv[1:dz] = np.exp(1.-1./sig[1:dz]**2) bpv[llmp1-1] = 0. apv[0] = 0. apv[1:dz+1] = pa*(sig[1:dz+1] - bpv[1:dz+1]) elif vert_sampling == 'read': # Read "ap" and "bp". First line is skipped (title line). "ap" # should be in Pa. First couple of values should correspond to # the surface, that is : "bp" should be in descending order. if not os.path.isfile('hybrid.txt'): print errormsg print ' ' + fname + ": parameters file 'hybrid.txt' does not exist!!" quit(-1) sfobj = ('hybrid.txt', 'r') # skip title line title = sfobj.readline() for l in range(dz+1): values = ncvar.readuce_space(sfobj.readline()) apv[l] = np.float(values[0]) bpv[l] = np.float(values[0]) sfobj.close() if apv[0] == 0. or apv[dz+1] == 0. or bpv[0] == 1. or bpv[dz+1] == 0.: print errormsg print ' ' + fname + ': bad ap or bp values !!' print ' k ap bp ___________' for k in range(dz+1): print k, apv[k], bpv[k] else: print errormsg print ' ' + fname + ': wrong value for vert_sampling:', vert_sampling quit(-1) nivsigs = np.arange(dz)*1.+1. nivsig = np.arange(llmp1)*1.+1. print ' ' + fname + ': k BP AP ________' for k in range(dz+1): print k, bpv[k], apv[k] print 'Niveaux de pressions approximatifs aux centres des' print 'couches calcules pour une pression de surface =', preff print 'et altitudes equivalentes pour une hauteur d echelle de ' print scaleheight,' km' for l in range(dz): dpres[l] = bpv[l] - bpv[l+1] pnivs[l] = 0.5*( apv[l]+bpv[l]*preff + apv[l+1]+bpv[l+1]*preff ) print ' PRESNIVS(', l, ')=', pnivs[l], ' Z ~ ', \ np.log(preff/pnivs[l])*scaleheight, ' DZ ~ ', \ scaleheight*np.log((apv[l]+bpv[l]*preff)/ \ np.max([apv[l+1]+bpv[l+1]*preff, 1.e-10])) print ' ' + fname + ': PRESNIVS [Pa]:', pnivs return pnivs, apv, bpv def presnivs_calc_noterre(hybrid, dz): """ From dyn3d/disvert_noterre.F calculation of vertical pressure levels hybrid= whether hydbrid coordinates have to be used dz= numbef of vertical layers """ fname = 'presnivs_calc_noterre' #----------------------------------------------------------------------- # .... Calculs de ap(l) et de bp(l) .... # ......................................... # # ..... pa et preff sont lus sur les fichiers start par dynetat0 ..... #----------------------------------------------------------------------- # llmp1 = dz + 1 pnivs = np.zeros((dz), dtype=np.float) s = np.zeros((dz), dtype=np.float) sig = np.zeros((dz+1), dtype=np.float) apv = np.zeros((dz+1), dtype=np.float) bpv = np.zeros((dz+1), dtype=np.float) # Ouverture possible de fichiers typiquement E.T. if os.path.isfile('easig.def'): #----------------------------------------------------------------------- # cas 1 on lit les options dans esasig.def: # ---------------------------------------- ofet = open('esasig.def', 'r') # Lecture de esasig.def : # Systeme peu souple, mais qui respecte en theorie # La conservation de l'energie (conversion Energie potentielle # <-> energie cinetique, d'apres la note de Frederic Hourdin... print '*****************************' print "WARNING reading 'esasig.def'" print '*****************************' for line in ofet: linevalues = ncvar.reduce_spaces(line) scaleheight = np.float(linevalues[0]) dz0 = np.float(linevalues[1]) dz1 = np.float(linevalues[2]) nhaut = np.float(linevalues[3]) ofet.close() dz0 = dz0/scaleheight dz1 = dz1/scaleheight sig1=(1.-dz1)/tanh(.5*(dz-1)/nhaut) esig=1. for l in range(19): esig=-np.log((1./sig1-1.)*np.exp(-dz0)/esig)/(dz-1.) csig=(1./sig1-1.)/(np.exp(esig)-1.) for l in range(1,dz): zz=csig*(np.exp(esig*(l-1.))-1.) sig[l] = 1./(1.+zz)*np.tanh(.5*(dz+1-l)/nhaut) sig[0] = 1. sig[dz] = 0. quoi = 1. + 2.* kappa s[dz-1] = 1. s[dz-2] = quoi if dz > 1: for ll in range(1, dz-2): l = dz+1 - ll quand = sig[l+1]/sig[l] s[l-1] = quoi*(1.-quand)*s[l] + quand*s[l+1] # snorm=(1.-.5*sig[1]+kappa*(1.-sig[1]))*s[0]+.5*sig[1]*s[1] for l in range(dz): s[l] = s[l]/ snorm elif os.path.isfile('z2sig.def'): fet = open('z2sig.def', 'r') #----------------------------------------------------------------------- # cas 2 on lit les options dans z2sig.def: # ---------------------------------------- print '****************************' print 'Reading z2sig.def' print '****************************' for line in ofet: linevalues = ncvar.reduce_spaces(line) scaleheight = np.float(linevalues[0]) for l in range(dz): zsig[l] = linevalues[l+1] ofet.close() sig[0] = 1. for l in range(1,dz): sig[l] = 0.5 * ( np.exp(-zsig[l]/scaleheight) + \ np.exp(-zsig[l-1]/scaleheight) ) sig[dz+1] = 0. #----------------------------------------------------------------------- else: print errormsg print ' ' + fname + ": didn't you forget something ???" print " We need file 'z2sig.def' ! (OR 'esasig.def')" quit(-1) #----------------------------------------------------------------------- nivsigs = np.arange(dz)*1. nivsig = np.arange(llmp1)*1. if hybrid: # use hybrid coordinates print "*********************************" print "Using hybrid vertical coordinates" print # Coordonnees hybrides avec mod for l in range(dz): newsig = sig_hybrid(sig[l],pa,preff) bpv[l] = np.exp(1.-1./(newsig**2)) apv[l] = pa * (newsig - bp[l] ) bp[llmp1-1] = 0. ap[llmp1-1] = 0. else: # use sigma coordinates print "********************************" print "Using sigma vertical coordinates" print # Pour ne pas passer en coordonnees hybrides for l in range(dz): apv[l] = 0. bpv[l] = sig[l] apv[llmp1-1] = 0. bp[llmp1-1] = 0. print 'BP________ ', bp print 'AP________ ', ap # Calcul au milieu des couches (llm = dz): # WARNING : le choix de placer le milieu des couches au niveau de # pression intermediaire est arbitraire et pourrait etre modifie. # Le calcul du niveau pour la derniere couche # (on met la meme distance (en log pression) entre P(llm) # et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est # Specifique. Ce choix est specifie ici ET dans exner_milieu.F for l in range(dz-1): aps[0:dz-1] = 0.5*( apv[0:dz-1]+apv[1:dz]) bps[0:dz-1] = 0.5*( bpv[0:dz-1]+bpv[1:dz]) if hybrid: aps[dz-1] = aps[dz-2]**2 / aps[dz-3] bps[dz-1] = 0.5*(bpv[dz-1] + bpv[dz]) else: bps[dz-1] = bps[dz-2]**2 / bps[dz-3] aps[dz-1] = 0. print 'BPs_______ ', bps print 'APs_______ ', aps for l in range(dz): psnivs[l] = aps[l]+bps[l]*preff psalt[l] = -scaleheight*np.log(presnivs[l]/preff) return psnivs, psalt def global_lonlat(dx,dy): """ Function to generate 2D matrices with the global longitude, latitudes dx, dy: dimension of the desired matrix >>> global_lonlat(5,5) array([[ 36., 108., 180., 252., 324.], [ 36., 108., 180., 252., 324.], [ 36., 108., 180., 252., 324.], [ 36., 108., 180., 252., 324.], [ 36., 108., 180., 252., 324.]]), array([[-72., -72., -72., -72., -72.], [-36., -36., -36., -36., -36.], [ 0., 0., 0., 0., 0.], [ 36., 36., 36., 36., 36.], [ 72., 72., 72., 72., 72.]])) """ fname = 'global_lonlat' longitude = np.zeros((dy,dx), dtype=np.float) latitude = np.zeros((dy,dx), dtype=np.float) for ix in range(dx): longitude[:,ix] = 360.*(1./2. + ix )/(dx) for iy in range(dy): latitude[iy,:] = 180.*(1./2. + iy )/(dy) - 90. return longitude, latitude def massdair(dx,dy,dz,p,airesurg): """c c ********************************************************************* c .... Calcule la masse d'air dans chaque maille .... c ********************************************************************* c c Auteurs : P. Le Van , Fr. Hourdin . c .......... c c .. p est un argum. d'entree pour le s-pg ... c .. masse est un argum.de sortie pour le s-pg ... c c .... p est defini aux interfaces des llm couches ..... c """ fname = 'massdair' masse = np.zeros((dz+1,dy+1,dx+1), dtype=np.float) # # # Methode pour calculer massebx et masseby . # ---------------------------------------- # # A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires # alpha1(i,j) calcule au point ( i+1/4,j-1/4 ) # alpha2(i,j) calcule au point ( i+1/4,j+1/4 ) # alpha3(i,j) calcule au point ( i-1/4,j+1/4 ) # alpha4(i,j) calcule au point ( i-1/4,j-1/4 ) # # Avec alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j) # # N.B . Pour plus de details, voir s-pg ... iniconst ... # # # # alpha4 . . alpha1 . alpha4 # (i,j) (i,j) (i+1,j) # # P . U . . P # (i,j) (i,j) (i+1,j) # # alpha3 . . alpha2 .alpha3 # (i,j) (i,j) (i+1,j) # # V . Z . . V # (i,j) # # alpha4 . . alpha1 .alpha4 # (i,j+1) (i,j+1) (i+1,j+1) # # P . U . . P # (i,j+1) (i+1,j+1) # # # # On a : # # massebx(i,j) = masse(i ,j) * ( alpha1(i ,j) + alpha2(i,j)) + # masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) ) # localise au point ... U (i,j) ... # # masseby(i,j) = masse(i,j ) * ( alpha2(i,j ) + alpha3(i,j ) + # masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1) # localise au point ... V (i,j) ... # # #======================================================================= for l in range (dz-1): for j in range(dy+1): for i in range(dx+1): masse[l,j,i] = airesurg[j,i] * ( p[l,j,i] - p[l+2,j,i] ) for j in range(1,dy): masse[l,j,dx] = masse[l,j-1,dx] return masse def geopot(dx, dy, dz, teta, pk, pks): """c======================================================================= c c Auteur: P. Le Van c ------- c c Objet: c ------ c c ******************************************************************* c .... calcul du geopotentiel aux milieux des couches ..... c ******************************************************************* c c .... l'integration se fait de bas en haut .... c c .. ngrid,teta,pk,pks,phis sont des argum. d'entree pour le s-pg .. c phi est un argum. de sortie pour le s-pg . c c======================================================================= """ fname = 'geopot' phis = np.zeros((dy+1,dx+1), dtype=np.float) phi = np.zeros((dz+1,dy+1,dx+1), dtype=np.float) print ' Lluis in ' + fname + ' shapes phis:',phis.shape,'phi:',phi.shape, \ 'teta:',teta.shape,'pks:',pks.shape,'pk:',pk.shape #----------------------------------------------------------------------- # calcul de phi au niveau 1 pres du sol ..... for j in range(dy): for i in range(dx): phi[0,j,i] = phis[j,i] + teta[0,j,i] * ( pks[j,i] - pk[0,j,i] ) # calcul de phi aux niveaux superieurs ....... for l in range(1,dz): for j in range(dy): for i in range(dx): phi[l,j,i] = phi[l-1,j,i] + 0.5 * ( teta[l,j,i]+teta[l-1,j,i] ) * \ ( pk[l-1,j,i]-pk[l,j,i] ) return phis, phi def dump2d(im,jm,nom_z): """ Function to create a dump 2d variable from dyn3d/dump2d.F """ fname = 'dump2d' z = np.zeros((im,jm), dtype=np.float) zmin = z[0,0] zllm = z[0,0] imin = 0 illm = 0 jmin = 0 jllm = 0 for j in range(jm): for i in range(im): if z[i,j] > zllm: illm=i jllm=j zllm=z[i,j] if z[i,j] < zmin: imin=i jmin=j zmin=z[i,j] print 'MIN:',zmin print 'MAX:',zllm if zllm > zmin: for j in range(jm): print int(10.*(z[:,j]-zmin)/(zllm-zmin)) return z def ugeostr(dx,dy,dz,phis,phi,rlatu,rlatv,cu): """! Calcul du vent covariant geostrophique a partir du champ de ! geopotentiel. ! We actually compute: (1 - cos^8 \phi) u_g ! to have a wind going smoothly to 0 at the equator. ! We assume that the surface pressure is uniform so that model ! levels are pressure levels. """ fname = 'ugeostr' ucov = np.zeros((dz,dy+1,dx+1), dtype=np.float) um = np.zeros((dz,dy), dtype=np.float) u = np.zeros((dz,dy,dx+1), dtype=np.float) print ' Lluis in ' + fname + ': shapes phis:',phis.shape,'phi:',phi.shape,'u:',u.shape for j in range(dy): if np.abs(np.sin(rlatv[j])) < 1.e-4: zlat = 1.e-4 else: zlat=rlatv[j] fact = np.cos(zlat) fact = fact*fact fact = fact*fact fact = fact*fact fact = (1.-fact)/ (2.*omeg*np.sin(zlat)*(rlatu[j+1]-rlatu[j])) fact = -fact/rad for l in range(dz): for i in range(dx): u[l,j,i] = fact*(phi[l,j+1,i]-phi[l,j,i]) um[l,j]=um[l,j]+u[l,j,i]/np.float(dx) um = dump2d(dz,dy,'Vent-u geostrophique') # calcul des champ de vent: for l in range(dz): for i in range(dx+1): ucov[l,0,i]=0. ucov[l,dy,i]=0. for j in range(1,dy): for i in range(dx): ucov[l,j,i] = 0.5*(u[l,j,i]+u[l,j-1,i])*cu[j,i] ucov[l,j,dx]=ucov[l,j,0] return ucov def RAN1(IDUM, Nvals): """ Function to generate Nvals random numbers from dyn3d/ran1.F IDUM= Random Seed Nvals= number of values """ fname = 'RAN1' R = np.zeros((Nvals), dtype=np.float) M1 = 259200 IA1 = 7141 IC1 = 54773 RM1 = 3.8580247E-6 M2 = 134456 IA2 = 8121 IC2 = 28411 RM2 = 7.4373773E-6 M3 = 243000 IA3 = 4561 IC3 = 51349 IFF = 0 if IDUM < 0 or IFF == 0: IFF = 1 IX1 = np.mod(IC1-IDUM,M1) IX1 = np.mod(IA1*IX1+IC1,M1) IX2 = np.mod(IX1,M2) IX1 = np.mod(IA1*IX1+IC1,M1) IX3 = np.mod(IX1,M3) for J in range(Nvals): IX1 = np.mod(IA1*IX1+IC1,M1) IX2 = np.mod(IA2*IX2+IC2,M2) R[J] = (np.float(IX1)+np.float(IX2)*RM2)*RM1 IDUM=1 IX1 = np.mod(IA1*IX1+IC1,M1) IX2 = np.mod(IA2*IX2+IC2,M2) IX3 = np.mod(IA3*IX3+IC3,M3) J = 1+(Nvals*IX3)/M3 if J > Nvals or J < 1: quit() ran1=R[J] R[J]=(np.float(IX1)+np.float(IX2)*RM2)*RM1 return ran1 def name_variables(filekind): """ Function to provide name of the variables and their atributes as function of the output type of file filekind= kind of file 'CF': CF-convention 'LMDZ': LMDZ style 'WRF': WRF style """ fname = 'name_variables' dimnames = {} varnames = {} # Standard dimensions dimn = ['x','y','z','t'] # Standard variables varn = ['lon', 'lat', 'lev', 'time', 'temp', 'tsfc', 'u10m', 'v10m', 'u', 'v', \ 'zsfc', 'geopot', 'psfc', 'pres', 'H2Ov', 'H2Ol'] # Standard variables' attribute names stdattrn = ['standard_name', 'long_name', 'units'] # Extra variables' attribute names # kextrattrn = [] # Dictionary with the values for each standard variable # kvarn['stdn'] = ['name', 'dims(list)', 'std_name', 'long_name', 'units', 'FillValue'] # None: No value kdimn = {} kvarn = {} kattrn = {} if filekind == 'CF': kdimn['x'] = 'x' kdimn['y'] = 'y' kdimn['z'] = 'z' kdimn['t'] = 'time' kvarn['lon'] = ['lon',['x'],'longitude','Longitude','degrees_east',None] kvarn['lat'] = ['lat',['y'],'latitude','Latitude','degrees_north',None] kvarn['lev'] = ['lev',['z'],'levels','Levels','-',None] kvarn['time'] = ['time',['t'],'time','Time', \ 'monutes since 1949-12-01 00:00:00', None] kvarn['temp'] = ['ta',['t','z','y','x'],'air_temperature','Air temperature', \ 'K',None] kvarn['tsfc'] = ['tas',['t','y','x'],'air_temperature','Air temperature','K',None] kvarn['u10m'] = ['uas',['t','y','x'],'eastward_wind','eastward wind','ms-1',None] kvarn['v10m'] = ['vas',['t','y','x'],'northward_wind','northward wind','ms-1',None] kvarn['u'] = ['ua',['t','z','y','x'],'eastward_wind','eastward wind','ms-1',None] kvarn['v'] = ['va',['t','z','y','x'],'northward_wind','northward wind','ms-1',None] kvarn['zsfc'] = ['zgs',['t','y','x'],'surface_geopotential_height', \ 'surface geopotential height','m2s-2',None] kvarn['geopot'] = ['zg',['t','z','y','x'],'geopotential_height','geopotential height', \ 'm2s-2',None] kvarn['psfc'] = ['ps',['t','z','y','x'],'surface_air_pressure','surface pressure','Pa',None] kvarn['pres'] = ['pres',['t','z','y','x'],'air_pressure','pressure','Pa',None] kvarn['H2Ov'] = ['r',['t','z','y','x'],'water_mixing_ratio','water mixing ratio','kgkg-1', \ None] kvarn['H2Ol'] = ['c',['t','z','y','x'],'condensed_water_mixing_ratio', \ 'condensed water mixing ratio','kgkg-1',None] kattrn['standard_name'] = 'standard_name' kattrn['long_name'] = 'long_name' kattrn['units'] = 'units' kextrattrn = [''] elif filekind == 'LMDZ': kdimn['x'] = 'x' kdimn['y'] = 'y' kdimn['z'] = 'presnivs' kdimn['t'] = 'time_counter' kvarn['lon'] = ['lon',['x'],'longitude','Longitude','degrees_east', None] kvarn['lat'] = ['lat',['y'],'latitude','Latitude','degrees_north', None] kvarn['lev'] = ['presnivs',['z'],'model_level_number','Vertical levels','Pa', \ None] kvarn['time'] = ['time_counter',['t'],'time','Time', \ 'seconds since 1980-01-01 00:24:00', None] kvarn['temp'] = ['temp',['t','z','y','x'],'Air temperature','Air temperature','K', \ 9.96921e+36] kvarn['tsfc'] = ['t2m',['t','y','x'],'Temperature 2m','Temperature 2m','K',9.96921e+36] kvarn['u10m'] = ['u10m',['t','y','x'],'Vent zonal 10m','Vent zonal 10m','m/s',9.96921e+36] kvarn['v10m'] = ['v10m',['t','y','x'],'Vent meridien 10m','Vent meridien 10m','m/s', \ 9.96921e+36] kvarn['u'] = ['vitu',['t','z','y','x'],'Zonal wind','Zonal wind','m/s',9.96921e+36] kvarn['v'] = ['vitv',['t','z','y','x'],'Meridional wind','Meridional wind','m/s',9.96921e+36] kvarn['zsfc'] = ['phis',['t','y','x'],'Surface geop.height','Surface geop.height', \ 'm2/s2',9.96921e+36] kvarn['geopot'] = ['geop',['t','z','y','x'],'Geopotential height','Geopotential height', \ 'm2/s2',9.96921e+36] kvarn['psfc'] = ['psol',['t','y','x'],'Surface Pressure','Surface Pressure','Pa', \ 9.96921e+36] kvarn['pres'] = ['pres',['t','z','y','x'],'Air pressure','Air pressure','Pa',9.96921e+36] kvarn['H2Ov'] = ['ovap',['t','z','y','x'],'Specific humidity','Specific humidity','kg/kg', \ 9.96921e+36] kvarn['H2Ol'] = ['oliq',['t','z','y','x'],'Condensed water','Condensed water','kg/kg', \ 9.96921e+36] kattrn['standard_name'] = 'standard_name' kattrn['long_name'] = 'long_name' kattrn['units'] = 'units' kextrattrn = ['coordinates'] elif filekind == 'WRF': kdimn['x'] = 'west_east' kdimn['y'] = 'south_north' kdimn['z'] = 'bottom_top' kdimn['t'] = 'Time' kvarn['lon'] = ['XLONG',['t','y','x'], None,'LONGITUDE, WEST IS NEGATIVE','degree_east', \ None] kvarn['lat'] = ['XLAT',['t','y','x'], None,'LATITUDE, SOUTH IS NEGATIVE','degree_north', \ None] kvarn['lev'] = ['ZNU',['t','z'],None,'eta values on half (mass) levels','',None] kvarn['time'] = ['Times',['t','DateStrLen'],None,None,None,None] kvarn['temp'] = ['T',['t','z','y','x'],None,'perturbation potential temperature (theta-t0)','K'\ ,None] kvarn['tsfc'] = ['T2',['t','y','x'],None,'TEMP at 2 M','K',None] kvarn['u10m'] = ['U10',['t','y','x'],None,'U at 10 M','m s-1',None] kvarn['v10m'] = ['V10',['t','y','x'],None,'V at 10 M','m s-1',None] kvarn['u'] = ['U',['t','z','y','x'],None,'x-wind component','m s-1',None] kvarn['v'] = ['V',['t','z','y','x'],None,'y-wind component','m s-1',None] kvarn['zsfc'] = [None,['t','y','x'],None,None,None,None] kvarn['geopot'] = ['PHB',['t','z','y','x'],None,'perturbation geopotential','m2 s-2',None] kvarn['psfc'] = [None,['t','y','x'],None,None,None,None] kvarn['pres'] = ['P',['t','z','y','x'],None,'perturbation pressure','Pa',None] kvarn['H2Ov'] = ['QVAPOR',['t','z','y','x'],None,'Water vapor mixing ratio','kg kg-1',None] kvarn['H2Ol'] = ['QCLOUD',['t','z','y','x'],None,'Cloud water mixing ratio','kg kg-1',None] kattrn['standard_name'] = None kattrn['long_name'] = 'description' kattrn['units'] = 'units' kextrattrn = ['FieldType','MemoryOrder','stagger'] # elif filekind == 'OTHER': # kdimn['x'] = '' # kdimn['y'] = '' # kdimn['z'] = '' # kdimn['t'] = '' # kvarn['lon'] = ['',['t','z','y','x'],'','','',] # kvarn['lat'] = ['',['t','z','y','x'],'','','',] # kvarn['lev'] = ['',['t','z','y','x'],'','','',] # kvarn['time'] = ['',['t','z','y','x'],'','','',] # kvarn['temp'] = ['',['t','z','y','x'],'','','',] # kvarn['tsfc'] = ['',['t','y','x'],'','','',] # kvarn['u10m'] = ['',['t','y','x'],'','','',] # kvarn['v10m'] = ['',['t','y','x'],'','','',] # kvarn['u'] = ['',['t','z','y','x'],'','','',] # kvarn['v'] = ['',['t','z','y','x'],'','','',] # kvarn['zsfc'] = ['',['t','y','x'],'','','',] # kvarn['geopot'] = ['',['t','z','y','x'],'','','',] # kvarn['psfc'] = ['',['t','y','x'],'','','',] # kvarn['pres'] = ['',['t','z','y','x'],'','','',] # kvarn['H2Ov'] = ['',['t','z','y','x'],'','','',] # kvarn['H2Ol'] = ['',['t','z','y','x'],'','','',] # kattrn['standard_name'] = '' # kattrn['long_name'] = '' # kattrn['units'] = '' # kextrattrn = [''] else: print errormsg print ' ' + fname + ": filekind '" + filekind + "' not ready !!" quit(-1) return kdimn, kvarn, kattrn, kextrattrn def generic_NetCDFfile(ncobj, dims, kfile, kdimns, kvarns, stdattrns, extrattrns): """ Function to fill a generic NetCDF file ncobj= NetCDF file object to which the variables have to be created kfile= kind of file 'CF': CF-convention 'LMDZ': LMDZ style 'WRF': WRF style dims= [dimx, dimy, dimz] dimensions of the file kdimns= dictionary for the specific names for the standard dimensions (x,y,z,t) kvarns= dictionary for the specific values for the standard variables kvarn['stdn'] = ['name', 'dims(list)', 'std_name', 'long_name', 'units', 'FillValue'] None: No value stdattrns= dictionary for the specific values for the variables' standard attributes extrattrns= list for the specific values for the variables' extra attributes """ fname = 'generic_NetCDFile' # Creation of dimensions ## newdim = ncobj.createDimension(kdimns['x'], dims[0]) newdim = ncobj.createDimension(kdimns['y'], dims[1]) newdim = ncobj.createDimension(kdimns['z'], dims[2]) newdim = ncobj.createDimension(kdimns['t'], None) if kfile == 'WRF': newdim = ncobj.createDimension('DateStrLen', 19) # Creation of variables ## for varn in kvarns.keys(): varvals = kvarns[varn] if varvals[0] is not None: dimsvar = [] for dimn in varvals[1]: if dimn == 'DateStrLen': dimsvar.append(dimn) else: dimsvar.append(kdimns[dimn]) if varvals[5] is not None: nerwvar = ncobj.createVariable(varvals[0], 'f4', tuple(dimsvar), \ fill_value=np.float(varvals[5])) else: newvar = ncobj.createVariable(varvals[0], 'f4', tuple(dimsvar)) # Attributes attrns = stdattrns.keys() for iattr in range(len(attrns)): attrn = attrns[iattr] if stdattrns[attrn] is not None and varvals[4+iattr] is not None: newvar.setncattr(stdattrns[attrn], varvals[4+iattr]) # Extra attributes for iattr in range(len(extrattrns)): attrn = extrattrns[iattr] if kfile == 'LMDZ': if attrn == 'coordinates': attrv = '' for din in varvals[1]: attrv = attrv + kdimns[din] + ' ' elif kfile == 'WRF': if attrn == 'FieldType': attrv = '104' elif attrn == 'MemoryOrder': attrv = '' Ndims = len(varvals[1]) for idim in range(Ndims-1,0,-1): if varvals[1][idim] != 't': attrv = attrv + varvals[1][idim].upper() elif attrn == 'stagger': staggeredvars = {} staggeredvars['U'] = 'X' staggeredvars['V'] = 'Y' staggeredvars['PH'] = 'Z' if ncvar.searchInlist(staggeredvars.keys(),varvals[0]): attrv = staggeredvars[varvals[0]] else: attrv = '' newvar.setncattr(extrattrns[iattr], attrv) ncobj.sync() return ####### ###### ##### #### ### ## # filekindsnames = "'" + ncvar.numVector_String(filekinds, "', '") + "'" parser = OptionParser() parser.add_option("-o", "--outputkind", type='choice', dest="okind", choices=filekinds, help="kind of output: " + filekindsnames, metavar="KIND") parser.add_option("-d", "--dimensions", dest="dims", help="dimensions of the initial conditions: dimx,dimy,dimz", metavar="VALUES") parser.add_option("-p", "--pressure_exner", dest="pexner", help="how as to b computed Exner pressure ('hybrid': hybrid coordinates, 'middle': middle layer)", metavar="VALUE") parser.add_option("-q", "--NWaterSpecies", dest="nqtot", help="Number of water species", metavar="VALUE") parser.add_option("-t", "--time", dest="time", help="time of the initial conditions ([YYYY][MM][DD][HH][MI][SS] format)", metavar="VALUE") parser.add_option("-z", "--z_levels", type='choice', dest="znivs", choices=['param', 'tropo', 'strato', 'read'], help="which kind of vertical levels have to be computed ('param', 'tropo', 'strato', 'read')", metavar="VAR") (opts, args) = parser.parse_args() ####### ####### ## MAIN ####### # dynamic variables # vcov, ucov: covariant winds # teta: potential temperature # q: advecting fields (humidity species) # ps: surface pressure # masse: air mass # phis: surface geopotential # Local: # p: pressure at the interface between layers (half-sigma) # pks: Exner function at the surface # pk: Exner functino at the half-sigma layers # pkf: Filtred Exner function at the half-sigma layers # phi: geopotential height # ddsin,zsig,tetapv,w_pv: auxiliar variables # tetastrat: potential temporeature in the stratosphere (K) # teta0,ttp,delt_y,delt_z,eps: constants for the T profile # k_f,k_c_a,k_c_s: calling constants # ok_geost: Initialisation geostrohic wind # ok_pv: Polar Vortex # phi_pv,dphi_pv,gam_pv: polar vortex constants dimx = int(opts.dims.split(',')[0]) dimy = int(opts.dims.split(',')[1]) dimz = int(opts.dims.split(',')[2]) nqtot = int(opts.nqtot) ofile = 'iniaqua.nc' print 'dimensions: ',dimx,', ',dimy,', ',dimz if opts.okind == 'CF': varnames = ['zg', 'ta', 'ua', 'va', 'pres', 'r'] timev = float(opts.time) # Reference time from 1949-12-01 00:00:00 UTC timev = ncvar.datetimeStr_conversion(opts.time,'YmdHMS','cfTime') dimnames = ['time', 'z', 'y', 'x'] elif opts.okind == 'WRF': varnames = ['Geopot', 'T', 'U', 'V', 'PRES', 'QVAPOR'] timev = ncvar.datetimeStr_conversion(opts.time,'YmdHMS','WRFdatetime') dimnames = ['Time', 'bottom_top', 'south_north', 'west_east'] # Constants ## llm = dimz ok_geost = True # Constants for Newtonian relaxation and friction k_f = 1. k_f = 1./(daysec*k_f) # cooling surface k_c_s=4. k_c_s=1./(daysec*k_c_s) # cooling free atm k_c_a=40. k_c_a=1./(daysec*k_c_a) # Constants for Teta equilibrium profile # mean Teta (S.H. 315K) teta0=315. # Tropopause temperature (S.H. 200K) ttp=200. # Deviation to N-S symmetry(~0-20K) eps=0. # Merid Temp. Gradient (S.H. 60K) delt_y=60. # Vertical Gradient (S.H. 10K) delt_z=10. # Polar vortex ok_pv = False # Latitude of edge of vortex phi_pv=-50. phi_pv=phi_pv*pi/180. # Width of the edge dphi_pv=5. dphi_pv=dphi_pv*pi/180. # -dT/dz vortex (in K/km) gam_pv=4. # For extra-terrestrial planets #presnivs, pseudoalt = presnivs_calc_noterre(opts.znivs, dimz) presnivs, ap, bp = presnivs_calc(opts.znivs, dimz) lon, lat = global_lonlat(dimy,dimx) lonu, latu = global_lonlat(dimy,dimx+1) lonv, latv = global_lonlat(dimy+1,dimx) # 2. Initialize fields towards which to relax ## knewt_t = np.zeros((dimz), dtype=np.float) kfrict = np.zeros((dimz), dtype=np.float) clat4 = np.zeros((dimy+1, dimx+1), dtype=np.float) # Friction knewt_g = k_c_a for l in range(dimz): zsig=presnivs[l]/preff knewt_t[l]=(k_c_s-k_c_a)*np.max([0.,(zsig-0.7)/0.3]) kfrict[l]=k_f*np.max([0.,(zsig-0.7)/0.3]) for j in 1,range(dimy+1): clat4[j,:]=np.cos(latu[j,0])**4 # Vertical profile tetajl = np.zeros((dimz, dimy+1, dimx), dtype=np.float) #theta = np.zeros((dimy+1, dimx+1), dtype=np.float) #thetarappel = np.zeros((dimz, dimy+1, dimx+1), dtype=np.float) theta = np.zeros((dimy, dimx), dtype=np.float) thetarappel = np.zeros((dimz, dimy, dimx), dtype=np.float) for l in range (dimz): zsig = presnivs[l]/preff tetastrat = ttp*zsig**(-kappa) tetapv = tetastrat if ok_pv and zsig < 0.1: tetapv = tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) ddsin = np.sin(latu) tetajl[l,:,:] = teta0-delt_y*ddsin+eps*ddsin-delt_z*(1.-ddsin*ddsin)*np.log(zsig) if planet_type == 'giant': tetajl[l,:,:] = teta0+(delt_y*((np.sin(latu*3.14159*eps+0.0001))**2) / \ ((latu*3.14159*eps+0.0001)**2))-delt_z*np.log(zsig) # Profile stratospheric isotherm (+vortex) for iy in range(dimy): for ix in range(dimx): w_pv=(1.-np.tanh((latu[iy,ix]-phi_pv)/dphi_pv))/2. tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv tetajl[l,iy,ix]=np.max([tetajl[l,iy,ix],tetastrat]) #for iz in range(dimz): # for iy in range(dimy+1): # tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,:] # # tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,dimx-1] thetarappel = tetajl.copy() # 3. Initialize fields (if necessary) # surface pressure if iflag_phys > 2: # specific value for CMIP5 aqua/terra planets # "Specify the initial dry mass to be equivalent to # a global mean surface pressure (101325 minus 245) Pa." press = np.ones((dimy+1, dimx+1), dtype=np.float)*101080. else: # use reference surface pressure press = np.ones((dimy+1, dimx+1), dtype=np.float)*preff # ground geopotential phiss = np.zeros((dimy+1, dimx+1), dtype=np.float) pres = pression(dimx, dimy, dimz, ap, bp, press) aire, apolnorth, apolsouth, airesurge, rlatitudeu, rlatitudev, cuwind, cvwind = \ inigeom(dimx, dimy) if opts.pexner == 'hybdrid': pks, pk, pkf, alpha, beta = exner_hyb(dimx, dimy, dimz, press, pres, aire, \ apolnorth, apolsouth) else: pks, pk, pkf = exner_milieu(dimx, dimy, dimz, press, pres, beta) masse = massdair(dimx,dimy,dimz,pres,airesurge) # bulk initialization of temperature theta = thetarappel.copy() # geopotential phisfc = np.zeros((dimy+1, dimx+1), dtype=np.float) phiall = np.zeros((dimz+1, dimy+1, dimx+1), dtype=np.float) phisfc, phiall = geopot(dimx,dimy,dimz,theta,pk,pks) # winds ucov = np.zeros((dimz, dimy, dimx), dtype=np.float) vcov = np.zeros((dimz, dimy, dimx), dtype=np.float) if ok_geost: ucov = ugeostr(dimx,dimy,dimz,phisfc,phiall,rlatitudeu,rlatitudev,cuwind) # bulk initialization of tracers q = np.zeros((dimz, dimy, dimx, nqtot), dtype=np.float) if planet_type == 'earth': # Earth: first two tracers will be water for i in range(nqtot): if i == 1: q[:,:,i] = 1.e-10 if i == 2: q[:,:,i] = 1.e-15 if i > 2: q[:,:,i] = 0. # add random perturbation to temperature idum = -1 zz = RAN1(idum,97) idum = 0 for l in range(dimz): for j in range(dimy): for i in range(dimx): theta[l,j,i] = theta[l,j,i]*(1.+0.005*RAN1(idum,97)) # maintain periodicity in longitude for l in range(dimz): for j in range(1,dimy): theta[l,j,dimx-1]=theta[l,j-1,dimx-1] ncf = NetCDFFile(ofile, 'w') # File structure creation filedimns, filevarns, filevarattr, filevarxtrattr = name_variables(opts.okind) generic_NetCDFfile(ncf,[dimx,dimy,dimz], opts.okind, filedimns, filevarns, \ filevarattr, filevarxtrattr) # File filling for varn in ncf.variables.keys(): varobj = ncf.variables[varn] if varn == filevarns['lon'][0]: varobj[:] = lon elif varn == filevarns['lat'][0]: varobj[:] = lat elif varn == filevarns['lev'][0]: varobj[:] = presnivs elif varn == filevarns['time'][0]: varobj[:] = timev elif varn == filevarns['temp'][0]: print 'Lluis shapes varobj:',varobj.shape,'theta:',theta.shape varobj[0,:,:,:] = theta[:] elif varn == filevarns['tsfc'][0]: varobj[:] = theta[0,:,:] elif varn == filevarns['u10m'][0]: varobj[:] = ucov[0,:,:] elif varn == filevarns['v10m'][0]: varobj[:] = vcov[0,:,:] elif varn == filevarns['u'][0]: varobj[0,:,:,:] = ucov[:] elif varn == filevarns['v'][0]: varobj[0,:,:,:] = vcov[:] elif varn == filevarns['zsfc'][0]: varobj[0,:,:] = phisfc[:] elif varn == filevarns['geopot'][0]: varobj[0,:,:] = phisall[:] elif varn == filevarns['psfc'][0]: varobj[0,:,:] = press[:] elif varn == filevarns['pres'][0]: varobj[0,:,:,:] = pres[:] elif varn == filevarns['H2Ov'][0]: varobj[0,:,:,:] = q[:,:,:0] elif varn == filevarns['H2Ol'][0]: varobj[0,:,:,:] = q[:,:,:1] ncf.sync() ncf.close() print main + ": successfull writing of file '" + ofile + "' !!"