source: branches/2015/dev_r5102_CNRS11_TestCase/NEMOGCM/TOOLS/TEST_CASES/CREATE_FILES/make_coordinates_gphi0.py @ 5346

#!/opt/local/bin/python
# python script to create coordinates files with e1=e2=constant, and gphit=gphiu=gphiv=gphif=0 (to do not have CORIOLIS)

import os,sys
from netCDF4 import Dataset as netcdf
import numpy as np
import math

res=0.25
resname='r025'

# output file
fcoord='coordinates_'+str(resname)+'_cst_gphi0.nc'

latmin=-30.
latmax=30.
latres=res

lonmin=0.
lonmax=180.
lonres=res

#LONGITUDE
nb_lonpoints=int((lonmax-lonmin)/lonres+1)
lontot=np.zeros(nb_lonpoints)

lon=lonmin
end=(lonmax-lonmin)

ii=0
while lon <= end :
   lontot[ii]=lon
   lon=lon+lonres
   ii=ii+1

# LATITUDE
# one point at equator needed
LAT1=int(latmax/latres+1)
LAT=2*LAT1-1
lat1=np.zeros(LAT1)
lat2=np.zeros(LAT1)
lat3=np.zeros(LAT)

lat=0.
ii=0
while lat <= latmax :
   lat1[ii]=lat
   lat=lat+latres
   ii=ii+1

lat=0.
ii=0
while lat >= latmin :
   lat2[ii]=lat
   lat=lat-latres
   ii=ii+1

ii=0
while ii < LAT1 :
   lat3[ii]=lat2[-1-ii]
   ii=ii+1

ii=0
while ii < LAT1 :
   lat3[LAT1-1+ii]=lat1[ii]
   ii=ii+1

print '\n    lat3:\n', lat3


print '   creating coordinates file  '+fcoord+'...'

LON=nb_lonpoints
nc=netcdf(fcoord,'w')
nc.createDimension('y',LAT)
nc.createDimension('x',LON)

cdfglamt=nc.createVariable('glamt','f',('y','x'))
cdfgphit=nc.createVariable('gphit','f',('y','x'))
cdfglamu=nc.createVariable('glamu','f',('y','x'))
cdfgphiu=nc.createVariable('gphiu','f',('y','x'))
cdfglamv=nc.createVariable('glamv','f',('y','x'))
cdfgphiv=nc.createVariable('gphiv','f',('y','x'))
cdfglamf=nc.createVariable('glamf','f',('y','x'))
cdfgphif=nc.createVariable('gphif','f',('y','x'))

cdfe1t=nc.createVariable('e1t','f',('y','x'))
cdfe1u=nc.createVariable('e1u','f',('y','x'))
cdfe1v=nc.createVariable('e1v','f',('y','x'))
cdfe1f=nc.createVariable('e1f','f',('y','x'))
cdfe2t=nc.createVariable('e2t','f',('y','x'))
cdfe2u=nc.createVariable('e2u','f',('y','x'))
cdfe2v=nc.createVariable('e2v','f',('y','x'))
cdfe2f=nc.createVariable('e2f','f',('y','x'))

# gphi = 0
for ii in range (LAT):
   cdfgphit[ii,:]=0.
   cdfgphiu[ii,:]=0.
   cdfgphiv[ii,:]=0.
   cdfgphif[ii,:]=0.


for ii in range (LON):  
   cdfglamt[:,ii]=lontot[ii]
   cdfglamu[:,ii]=lontot[ii] + lonres/2
   cdfglamv[:,ii]=lontot[ii]
   cdfglamf[:,ii]=lontot[ii] + lonres/2

# from DOM/phycst.F90
ra=6371229.                     # earth radius (meter)
rpi=math.pi
rad=rpi / 180.     # conversion from degre into radian

# e2 = constant
cdfe2t[:,:]=latres*rad*ra
cdfe2u[:,:]=latres*rad*ra
cdfe2v[:,:]=latres*rad*ra
cdfe2f[:,:]=latres*rad*ra

# e1 = constant
cdfe1t[:,:]=lonres*rad*ra
cdfe1u[:,:]=lonres*rad*ra
cdfe1v[:,:]=lonres*rad*ra
cdfe1f[:,:]=lonres*rad*ra


nc.close()

#sys.exit()