global_ocean.90x40x15/input_ad.kapredi/get_grdchk_location.py

61b91e10c1 Ivan*0001 #!/usr/bin/env python
                0002 # -*- coding: iso-8859-15 -*-
                0003 ######################## -*- coding: utf-8 -*-
                0004 """Usage: get_grdchk_location.py
                0005 purpose: find location of grdchk points from "data.grdchk" settings
                0006 """
                0007
                0008 #
                0009 # this is the python version of the matlab script get_grdchk_location.py
                0010 #
                0011 import numpy as np
                0012 import sys
                0013 sys.path.append('../../../utils/python/MITgcmutils')
                0014 import MITgcmutils as mit
                0015
                0016 # from a model run
                0017 hFacC=mit.rdmds('hFacC')
                0018 nr,ny,nx = hFacC.shape
                0019
                0020 # from SIZE.h :
                0021 sNx, sNy = nx//2, ny//2
                0022
                0023 # from data.grdchk :
                0024 iGloPos  = 31
                0025 jGloPos  = 7
                0026 kGloPos  = 1
                0027 iGloTile = 2
                0028 jGloTile = 2
                0029
                0030 # from data.grdchk :
                0031 # some parameters to test:
                0032 nsteps = [450, 560, 553]
                0033 nend = 1700
                0034
                0035 # count number of wet-points in this tile(iGloTile,jGloTile):
                0036 nSx, nSy = nx//sNx, ny//sNy
                0037
                0038 msk=np.ceil(hFacC).reshape((nr,nSy,sNy,nSx,sNx))
                0039 mskTile=np.squeeze(msk[:,iGloTile-1,:,jGloTile-1,:])
                0040
                0041 cumMsk=np.cumsum(mskTile.reshape((sNx*sNy*nr,)))
                0042
                0043 # find location:
                0044 ic1,jc1,kc1 = iGloPos,jGloPos,kGloPos
                0045
                0046 ijk1=ic1+(jc1-1)*sNx+(kc1-1)*sNx*sNy
                0047 yy1=cumMsk[ijk1-1]
                0048
                0049 for nstep in nsteps:
                0050     #-- find the location of all grdchk points:
                0051     print('In tile (%i,%i), nstep = %4i , (ic,jc,kc) = ( %2i , %2i , %2i ) :'%(
                0052         iGloTile,jGloTile, nstep,ic1,jc1,kc1))
                0053     for n in np.arange(0,nend+1,nstep):
                0054         I=np.where(cumMsk==(yy1+n))[0]
                0055         ii=I[0];               kg=1+np.floor(ii/(sNx*sNy))
                0056         ii=np.mod(ii,sNx*sNy); jg=1+np.floor(ii/sNx)
                0057         ii=np.mod(ii,sNx);     ig=1+ii
                0058         print(' #%6i : i,j,k = %2i , %2i , %2i'%(n,ig,jg,kg))
                0059
                0060     print(' ------------')