Changeset 215 in lmdz_wrf

Timestamp:

Jan 13, 2015, 1:11:32 PM (10 years ago)

Author:

lfita

Message:

Working but not tested version of the script!

File:

• trunk/tools/iniaqua.py (modified) (27 diffs)

trunk/tools/iniaqua.py

@@ -19 +19 @@
 filekinds = ['CF', 'WRF']
 
-## g.e. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z tropo
+## g.e. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z tropo -q 2
 def fxy(dx, dy): 
     """! 
@@ -1395 +1395 @@
     print rlatuu[0:dy+2]
 
-    return aire
+    return aire, apoln, apols, airesurg, rlatu, rlatv, cu, cv
 
 def SSUM(n,sx,incx):
@@ -1415 +1415 @@
     return ssumv
 
-def SCOPY(n,sx,incx,sy,incy):
+def SCOPY(ddx,ddy,ddz,sx,incx,sy,incy):
     """ Obsolete function to copy matrix values
       from dyn3d/cray.F
     """
-    iy = 1
-    ix = 1
-
-    for i in range(n):
-        sy[iy] = sx[ix]
-        ix = ix+incx
-        iy = iy+incy
+    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):
+def exner_hyb (dx, dy, dz, psv, pv, aire, apoln, apols):
     """c
        c     Auteurs :  P.Le Van  , Fr. Hourdin  .
@@ -1459 +1469 @@
     fname = 'exner_hyb'
 
-    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((dy, dx), dtype=np.float)
-    pps = np.zeros((dy, dx), dtype=np.float)
-
-    ip1jm = (dx+1)*dy
-
-    if dz == 1:       
+    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
@@ -1485 +1496 @@
     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):
-            pks[j,i] = cpp * ( ps[j,i]/preff ) ** kappa
-
-    for j in range(dy):
-        for i in range(dx):
-            ppn[j,i] = aire[j,i] * pks[j,i]
-            pps[j,i] = aire[j,i+ip1jm] * pks[ij+ip1jm]
-
-    xpn = SSUM(iim,ppn,1) /apoln
-    xps = SSUM(iim,pps,1) /apols
-
-    for ij in range(iip1):
-        pks[ij] = xpn
-        pks[ij+ip1jm] = xps
+            pksv[j,i] = xpn[i]
+            pksv[dy-1,i] = xps[i]
 #
 #
 #    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
 #
-    for ij in range(ngrid):
-        alpha[ij,llm] = 0.
-        beta [ij,llm] = 1./ unpl2k
+    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 (llm-1,1,-1):
-
-        for ij in range(ngrid):
-            dellta = p[ij,l]* unpl2k + p[ij,l+1]* ( beta[ij,l+1]-unpl2k )
-            alpha[ij,l] = -p[ij,l+1] / dellta * alpha[ij,l+1]
-            beta[ij,l] = p[ij,l] / dellta   
+    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 ij in range(ngrid):
-        pk[ij,0] = ( p[ij,0]*pks[ij] - 0.5*alpha[ij,1]*p[ij,1] )                     \
-          *(  p[ij,0]* (1.+kappa) + 0.5*( beta[ij,1]-unpl2k )* p[ij,1] )
+    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(1,llm):
-        for ij in range(ngrid):
-            pk[ij,l] = alpha[ij,l] + beta[ij,l] * pk[ij,l-1]
-#
-#
-    pkfv = SCOPY ( ngrid * llm, pk, 1, pkfv, 1 )
+    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, aplhav, betav
+    return pksv, pkv, pkfv, alphav, betav
 
 def exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ):
@@ -1599 +1614 @@
 #         Compute pks(:),pk(:),pkf(:)
         
-        for ij in range(ngrid):
-            pks[ij] = (cpp/preff) * ps[ij]
-            pk[ij,1] = .5*pks[ij]
+        for j,i in range(ngrid):
+            pks[j,i] = (cpp/preff) * ps[j,i]
+            pk[j,i,1] = .5*pks[j,i]
           
-        pkf = SCOPY( ngrid * llm, pk, 1, pkf, 1 )
+        pkf = SCOPY(dx,dy,dz, pk, 1, pkf, 1 )
 # We do not filter for iniaqua
 #        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
@@ -1616 +1631 @@
 #     -------------
    
-    for ij in range(ngrid):
-        pks[ij] = cpp * ( ps[ij]/preff ) ** kappa
-
-    for ij in range(iim):
-        ppn[ij] = aire[ij] * pks[ij]
-        pps[ij] = aire[ij+ip1jm] * pks[ij+ip1jm]
+    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 ij in range(iip1):
-        pks[ij] = xpn
-        pks[ij+ip1jm] = xps
+    for j,i in range(iip1):
+        pks[j,i] = xpn
+        pks[j,i+ip1jm] = xps
 
 #
@@ -1637 +1652 @@
     dum1 = cpp * (2*preff)**(-kappa) 
     for l in range(llm-1):
-        for ij in range(ngrid):
-            pk[ij,l] = dum1 * (p[ij,l] + p[ij,l+1])**kappa
+        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 ..
@@ -1644 +1659 @@
 #    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
 
-    for ij in range(ngrid):
-        pk[ij,llm] = pk[ij,llm-1]**2 / pk[ij,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( ngrid * llm, pk, 1, pkf, 1 )
+    pkf = SCOPY( dx,dy,dz, pk, 1, pkf, 1 )
 
 # We do not filter iniaqua
@@ -1657 +1672 @@
 #    --------------------------------
     for l in range(1, llm):
-        for ij in range(ngrid):
-            beta[ij,l] = pk[ij,l] / pk[ij,l-1]
+        for j,i in range(ngrid):
+            beta[j,i,l] = pk[j,i,l] / pk[j,i,l-1]
 
     return pksv, pkv, pkfv
@@ -1676 +1691 @@
     fname = 'pression'
       
-    press = np.zeros((dz+1, dy, dx), dtype=np.float)
-
-    print 'shape psv:',psv.shape,'press:',press.shape,'ap:',apv.shape,'bp:',bpv.shape
+    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
@@ -2099 +2113 @@
     return longitude, latitude
 
-def massdair(p):
+def massdair(dx,dy,dz,p,airesurg):
     """c
        c *********************************************************************
@@ -2116 +2130 @@
     fname = 'massdair'
 
-    masse = np.zeros((ip1jmp1,llm), dtype=np.float)
+    masse = np.zeros((dz+1,dy+1,dx+1), dtype=np.float)
 #
 #
@@ -2167 +2181 @@
 #=======================================================================
 
-    for l in range (llm):
-        for ij in range(ip1jmp1):
-            masse[ij,l] = airesurg[ij] * ( p[ij,l] - p[ij,l+1] )
-
-        for ij in range(ip1jmp1,iip1):
-            masse[ij+iim,l] = masse[ij,l]
+    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(ngrid, teta, pk, pks):
+def geopot(dx, dy, dz, teta, pk, pks):
     """c=======================================================================
        c
@@ -2198 +2213 @@
     fname = 'geopot'
 
-    phis = np.zeros((ngrid), dtype=np.float)
-    phi = np.zeros((ngrid,llm), dtype=np.float)
+    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 ij in range(ngrid):
-        phi[ij,1] = phis[ij] + teta[ij,0] * ( pks[ij] - pk[ij,0] )
+    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,llm):
-        for ij in range(ngrid):
-            phi[ij,l] = phi[ij,l-1] + 0.5 * ( teta[ij,l]+teta[ij,l-1] ) *            \
-              ( pk[ij,l-1]-pk[ij,l] )
+    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 
@@ -2252 +2272 @@
     return z
 
-def ugeostr(phi):
+def ugeostr(dx,dy,dz,phis,phi,rlatu,rlatv,cu):
     """! Calcul du vent covariant geostrophique a partir du champ de
        ! geopotentiel.
@@ -2261 +2281 @@
     """
     fname = 'ugeostr'
-    ucov = np.zeros((iip1,jjp1,llm), dtype=np.float)
-    um = np.zeros((jjm,llm), dtype=np.float)
-    u = np.zeros((iip1,jjm,llm), dtype=np.float)
-
-    for j in range(jjm):
+    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)
+            zlat=rlatv[j]
 
         fact = np.cos(zlat)
@@ -2275 +2297 @@
         fact = fact*fact
         fact = fact*fact
-        fact = (1.-fact)/ (2.*omeg*sin(zlat)*(rlatu[j+1]-rlatu[j]))
+        fact = (1.-fact)/ (2.*omeg*np.sin(zlat)*(rlatu[j+1]-rlatu[j]))
         fact = -fact/rad
-        for l in range(llm):
-            for i in range(iim):
-                u[i,j,l] = fact*(phi[i,j+1,l]-phi[i,j,l])
-                um[j,l]=um[j,l]+u[i,j,l]/np.float(iim)
-
-    um = dump2d(jjm,llm,'Vent-u geostrophique')
+        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(llm):
-        for i in range(iip1):
-            ucov[i,1,l]=0.
-            ucov[i,jjp1,l]=0.
-        for j in range(1,jjm):
-            for i in range(iim):
-                ucov[i,j,l] = 0.5*(u[i,j,l]+u[i,j-1,l])*cu[i,j]
-
-        ucov[iip1,j,l]=ucov[0,j,l]
+    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
@@ -2305 +2327 @@
     """
     fname = 'RAN1'
-    Nvals = 97
 
     R = np.zeros((Nvals), dtype=np.float)
@@ -2320 +2341 @@
     IA3 = 4561
     IC3 = 51349
+    IFF = 0
 
     if IDUM < 0 or IFF == 0:
@@ -2358 +2380 @@
   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")
@@ -2396 +2420 @@
 dimy = int(opts.dims.split(',')[1])
 dimz = int(opts.dims.split(',')[2])
+nqtot = int(opts.nqtot)
 
 ofile = 'iniaqua.nc'
@@ -2473 +2498 @@
 # Vertical profile
 tetajl =  np.zeros((dimz, dimy+1, dimx), dtype=np.float)
-teta = np.zeros((dimy+1, dimx+1), dtype=np.float)
-tetarappel = np.zeros((dimz, dimy+1, dimx+1), dtype=np.float)
+theta = np.zeros((dimy+1, dimx+1), dtype=np.float)
+thetarappel = np.zeros((dimz, dimy+1, dimx+1), dtype=np.float)
 
 for l in range (dimz):
@@ -2501 +2526 @@
 #
 #    tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,dimx-1]
-tetarappel = tetajl
+thetarappel = tetajl.copy()
 
 # 3. Initialize fields (if necessary)
@@ -2510 +2535 @@
 # "Specify the initial dry mass to be equivalent to
 #  a global mean surface pressure (101325 minus 245) Pa."
-    ps = np.ones((dimy, dimx), dtype=np.float)*101080.  
+    press = np.ones((dimy+1, dimx+1), dtype=np.float)*101080.  
 else:
 # use reference surface pressure
-    ps = np.ones((dimy, dimx), dtype=np.float)*preff
+    press = np.ones((dimy+1, dimx+1), dtype=np.float)*preff
         
 # ground geopotential
-phiss =  np.zeros((dimy, dimx), dtype=np.float)
-
-p = pression(dimx, dimy, dimz, ap, bp, ps)
-
-aire = inigeom(dimx, dimy)
+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, ps, p, aire)
+    pks, pk, pkf, alpha, beta = exner_hyb(dimx, dimy, dimz, press, pres, aire,       \
+      apolnorth, apolsouth)
 else:
-    pks, pk, pkf = exner_milieu(ip1jmp1, dimx, dimy, dimz, ps, p, beta)
-
-masse = massdair(p)
+    pks, pk, pkf = exner_milieu(ip1jmp1, dimx, dimy, dimz, press, pres, beta)
+
+masse = massdair(dimx,dimy,dimz,pres,airesurge)
 
 # bulk initialization of temperature
-teta = tetarappel.copy()
+theta = thetarappel.copy()
 
 # geopotential
-phis =  np.zeros((dimy, dimx), dtype=np.float)
-phi =  np.zeros((dimz, dimy, dimx), dtype=np.float)
-
-phis, phi = geopot(ip1jmp1,teta,pk,pks)
+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
@@ -2542 +2569 @@
 vcov =  np.zeros((dimz, dimy, dimx), dtype=np.float)
 
-
 if ok_geost:
-    ucov = ugeostr(phi)
+    ucov = ugeostr(dimx,dimy,dimz,phisfc,phiall,rlatitudeu,rlatitudev,cuwind)
 
 # bulk initialization of tracers
+print 'Lluis nqtot:',nqtot
 q = np.zeros((dimz, dimy, dimx, nqtot), dtype=np.float)
 
@@ -2558 +2585 @@
 # add random perturbation to temperature
 idum  = -1
-zz = RAN1(idum)
+zz = RAN1(idum,97)
 idum  = 0
-for l in range(llm):
-    for ij in range(iip2,ip1jm):
-        teta[ij,l] = teta[ij,l]*(1.+0.005*RAN1(idum))
+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(llm):
-    for ij in range(0,ip1jmp1,iip1):
-        teta[ij+iim,l]=teta[ij,l]
+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')