![(please configure the [header_logo] section in trac.ini)](/LMDZ_WRF/chrome/site/your_project_logo.png){kind=logo}
Changeset 2311 in lmdz_wrf
- Timestamp:
- Feb 1, 2019, 7:59:22 PM (6 years ago)
- Location:
- trunk/tools
- Files:
-
- 2 edited
-
module_scientific.f90 (modified) (2 diffs)
-
nc_var_tools.py (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/tools/module_scientific.f90
r2308 r2311 5701 5701 DO s3 =1, ds3 5702 5702 Nin = Ngridsin(s1,s2,s3) 5703 IF (ALLOCATED(gin)) DEALLOCATE(gin)5704 ALLOCATE(gin(Nin,2))5705 IF (ALLOCATED(pin)) DEALLOCATE(pin)5706 ALLOCATE(pin(Nin))5707 IF (ALLOCATED(vin)) DEALLOCATE(vin)5708 ALLOCATE(vin(Nin,di3))5709 IF (ALLOCATED(svin)) DEALLOCATE(svin)5710 ALLOCATE(svin(Nin))5711 gin = gridsin(s1,s2,s3,1:Nin,:)5712 pin = percentages(s1,s2,s3,1:Nin)5713 5714 ! Getting the values5715 DO iv=1, Nin5716 i1 = gin(iv,1)5717 i2 = gin(iv,2)5718 vin(iv,:) = varin(i1,i2,:)5719 END DO5720 5721 5703 ! Computing along d3 5722 DO i3=1, di3 5723 minv = fillVal64 5724 maxv = -fillVal64 5725 meanv = zeroRK 5726 mean2v = zeroRK 5727 stdv = zeroRK 5728 minv = MINVAL(vin(:,i3)) 5729 maxv = MAXVAL(vin(:,i3)) 5730 meanv = SUM(vin(:,i3)*pin) 5731 mean2v = SUM(vin(:,i3)**2*pin) 5732 DO iv=1,Nin 5733 stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2 5704 IF (Nin > 1) THEN 5705 IF (ALLOCATED(gin)) DEALLOCATE(gin) 5706 ALLOCATE(gin(Nin,2)) 5707 IF (ALLOCATED(pin)) DEALLOCATE(pin) 5708 ALLOCATE(pin(Nin)) 5709 IF (ALLOCATED(vin)) DEALLOCATE(vin) 5710 ALLOCATE(vin(Nin,di3)) 5711 IF (ALLOCATED(svin)) DEALLOCATE(svin) 5712 ALLOCATE(svin(Nin)) 5713 gin = gridsin(s1,s2,s3,1:Nin,:) 5714 pin = percentages(s1,s2,s3,1:Nin) 5715 5716 ! Getting the values 5717 DO iv=1, Nin 5718 i1 = gin(iv,1) 5719 i2 = gin(iv,2) 5720 vin(iv,:) = varin(i1,i2,:) 5734 5721 END DO 5735 stdv = SQRT(stdv) 5736 svin = vin(:,i3) 5737 CALL SortR_K(svin, Nin) 5738 medv = svin(INT(Nin/2)) 5739 varout(s1,s2,s3,i3,1) = minv 5740 varout(s1,s2,s3,i3,2) = maxv 5741 varout(s1,s2,s3,i3,3) = meanv 5742 varout(s1,s2,s3,i3,4) = mean2v 5743 varout(s1,s2,s3,i3,5) = stdv 5744 varout(s1,s2,s3,i3,6) = medv 5745 varout(s1,s2,s3,i3,7) = Nin*1. 5746 END DO 5722 DO i3=1, di3 5723 minv = fillVal64 5724 maxv = -fillVal64 5725 meanv = zeroRK 5726 mean2v = zeroRK 5727 stdv = zeroRK 5728 minv = MINVAL(vin(:,i3)) 5729 maxv = MAXVAL(vin(:,i3)) 5730 meanv = SUM(vin(:,i3)*pin) 5731 mean2v = SUM(vin(:,i3)**2*pin) 5732 DO iv=1,Nin 5733 stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2 5734 END DO 5735 stdv = SQRT(stdv) 5736 svin = vin(:,i3) 5737 CALL SortR_K(svin, Nin) 5738 medv = svin(INT(Nin/2)) 5739 varout(s1,s2,s3,i3,1) = minv 5740 varout(s1,s2,s3,i3,2) = maxv 5741 varout(s1,s2,s3,i3,3) = meanv 5742 varout(s1,s2,s3,i3,4) = mean2v 5743 varout(s1,s2,s3,i3,5) = stdv 5744 varout(s1,s2,s3,i3,6) = medv 5745 varout(s1,s2,s3,i3,7) = Nin*1. 5746 END DO 5747 ELSE 5748 varout(s1,s2,s3,:,1) = varin(i1,i2,:) 5749 varout(s1,s2,s3,:,2) = varin(i1,i2,:) 5750 varout(s1,s2,s3,:,3) = varin(i1,i2,:) 5751 varout(s1,s2,s3,:,4) = varin(i1,i2,:)*varin(i1,i2,:) 5752 varout(s1,s2,s3,:,5) = zeroRK 5753 varout(s1,s2,s3,:,6) = varin(i1,i2,:) 5754 varout(s1,s2,s3,:,7) = Nin*1. 5755 END IF 5747 5756 END DO 5748 5757 END DO … … 5817 5826 DO s4 =1, ds4 5818 5827 Nin = Ngridsin(s1,s2,s3,s4) 5819 IF (ALLOCATED(gin)) DEALLOCATE(gin) 5820 ALLOCATE(gin(Nin,2)) 5821 IF (ALLOCATED(pin)) DEALLOCATE(pin) 5822 ALLOCATE(pin(Nin)) 5823 IF (ALLOCATED(vin)) DEALLOCATE(vin) 5824 ALLOCATE(vin(Nin,di3)) 5825 IF (ALLOCATED(svin)) DEALLOCATE(svin) 5826 ALLOCATE(svin(Nin)) 5827 gin = gridsin(s1,s2,s3,s4,1:Nin,:) 5828 pin = percentages(s1,s2,s3,s4,1:Nin) 5829 5830 minv = fillVal64 5831 maxv = -fillVal64 5832 meanv = zeroRK 5833 mean2v = zeroRK 5834 stdv = zeroRK 5835 5836 ! Getting the values 5837 DO iv=1, Nin 5838 i1 = gin(iv,1) 5839 i2 = gin(iv,2) 5840 vin(iv,:) = varin(i1,i2,:) 5841 END DO 5842 5843 ! Computing along d3 5844 DO i3=1, di3 5845 minv = fillVal64 5846 maxv = -fillVal64 5847 meanv = zeroRK 5848 mean2v = zeroRK 5849 stdv = zeroRK 5850 5851 minv = MINVAL(vin(:,i3)) 5852 maxv = MAXVAL(vin(:,i3)) 5853 meanv = SUM(vin(:,i3)*pin) 5854 mean2v = SUM(vin(:,i3)**2*pin) 5855 DO iv=1,Nin 5856 stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2 5828 IF (Nin > 1) THEN 5829 IF (ALLOCATED(gin)) DEALLOCATE(gin) 5830 ALLOCATE(gin(Nin,2)) 5831 IF (ALLOCATED(pin)) DEALLOCATE(pin) 5832 ALLOCATE(pin(Nin)) 5833 IF (ALLOCATED(vin)) DEALLOCATE(vin) 5834 ALLOCATE(vin(Nin,di3)) 5835 IF (ALLOCATED(svin)) DEALLOCATE(svin) 5836 ALLOCATE(svin(Nin)) 5837 gin = gridsin(s1,s2,s3,s4,1:Nin,:) 5838 pin = percentages(s1,s2,s3,s4,1:Nin) 5839 5840 ! Getting the values 5841 DO iv=1, Nin 5842 i1 = gin(iv,1) 5843 i2 = gin(iv,2) 5844 vin(iv,:) = varin(i1,i2,:) 5857 5845 END DO 5858 stdv = SQRT(stdv) 5859 svin = vin(:,i3) 5860 CALL SortR_K(svin, Nin) 5861 medv = svin(INT(Nin/2)) 5862 varout(s1,s2,s3,s4,i3,1) = minv 5863 varout(s1,s2,s3,s4,i3,2) = maxv 5864 varout(s1,s2,s3,s4,i3,3) = meanv 5865 varout(s1,s2,s3,s4,i3,4) = mean2v 5866 varout(s1,s2,s3,s4,i3,5) = stdv 5867 varout(s1,s2,s3,s4,i3,6) = medv 5868 varout(s1,s2,s3,s4,i3,7) = Nin*1. 5869 END DO 5846 ! Computing along d3 5847 DO i3=1, di3 5848 minv = fillVal64 5849 maxv = -fillVal64 5850 meanv = zeroRK 5851 mean2v = zeroRK 5852 stdv = zeroRK 5853 5854 minv = MINVAL(vin(:,i3)) 5855 maxv = MAXVAL(vin(:,i3)) 5856 meanv = SUM(vin(:,i3)*pin) 5857 mean2v = SUM(vin(:,i3)**2*pin) 5858 DO iv=1,Nin 5859 stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2 5860 END DO 5861 stdv = SQRT(stdv) 5862 svin = vin(:,i3) 5863 CALL SortR_K(svin, Nin) 5864 medv = svin(INT(Nin/2)) 5865 varout(s1,s2,s3,s4,i3,1) = minv 5866 varout(s1,s2,s3,s4,i3,2) = maxv 5867 varout(s1,s2,s3,s4,i3,3) = meanv 5868 varout(s1,s2,s3,s4,i3,4) = mean2v 5869 varout(s1,s2,s3,s4,i3,5) = stdv 5870 varout(s1,s2,s3,s4,i3,6) = medv 5871 varout(s1,s2,s3,s4,i3,7) = Nin*1. 5872 END DO 5873 ELSE 5874 varout(s1,s2,s3,s4,:,1) = varin(i1,i2,:) 5875 varout(s1,s2,s3,s4,:,2) = varin(i1,i2,:) 5876 varout(s1,s2,s3,s4,:,3) = varin(i1,i2,:) 5877 varout(s1,s2,s3,s4,:,4) = varin(i1,i2,:)*varin(i1,i2,:) 5878 varout(s1,s2,s3,s4,:,5) = zeroRK 5879 varout(s1,s2,s3,s4,:,6) = varin(i1,i2,:) 5880 varout(s1,s2,s3,s4,:,7) = Nin*1. 5881 END IF 5870 5882 END DO 5871 5883 END DO -
trunk/tools/nc_var_tools.py
r2308 r2311 28488 28488 aa = aa + areas[iy,ix] 28489 28489 pp = pp + percens[iv,j,i] 28490 if np.isnan(slicearea): 28491 print errormsg 28492 print ' ' + fname + ': Nan slice area for slice:', \ 28493 j, i, 'N grids:', Ngridsin[j,i],' !!' 28494 print ' #grid iy,ix grid_area grid_percen______' 28495 for iv in range(Ngridsin[j,i]): 28496 ix = gridsin[1,iv,j,i] 28497 iy = gridsin[0,iv,j,i] 28498 print iv, iy, ix, areas[iy,ix], percens[j,i] 28499 quit(-1) 28490 28500 aanewvar[j,i]= slicearea 28491 28501 # This should be true for the regular lon/lat projections … … 28623 28633 anewvar[iv,jB,iB,jA,iA]= (aA*pA)*(aB*pB) 28624 28634 slicearea = slicearea + (aA*pA)*(aB*pB) 28635 if np.isnan(slicearea): 28636 print errormsg 28637 print ' ' + fname + ': Nan slice area for slice:', \ 28638 jA, iA, jB, iB, 'N grids:', Nin,' !!' 28639 print ' #grid A grid_area A grid_percen B grid_area ' + \ 28640 'B grid_percen______' 28641 for iv in range(Nin): 28642 print iv, aA, pA, aB, pB 28643 quit(-1) 28625 28644 aanewvar[jB,iB,jA,iA]= slicearea 28626 28645 … … 28786 28805 ' slice:', jA, iA, jB, iB, 'N grids:', Nin,' !!' 28787 28806 print ' slice total area:', aanewvar[jA,iA,jB,iB] 28788 print ' #grid grid_area grid_percen _______'28807 print ' #grid grid_area slice_area grid_percen ______' 28789 28808 for iv in range(Nin): 28790 print iv, anewvar[iv,jA,iA,jB,iB], sgpa[iv,jA,iA,jB,iB] 28809 print iv, anewvar[iv,jA,iA,jB,iB], \ 28810 aanewvar[jA,iA,jB,iB], sgpa[iv,jA,iA,jB,iB] 28791 28811 quit(-1) 28792 28812 onewnc.sync() … … 28833 28853 print ' ' + varn + ' ...' 28834 28854 ovar = onc.variables[varn] 28835 varv = ovar[:]28836 28855 vdimns = list(ovar.dimensions) 28837 28856 vshape = list(ovar.shape) … … 28903 28922 voutt = None 28904 28923 varvt = None 28905 varvt = varv.transpose() 28906 if len(vshape) == 3 and len(sN.shape) == 4: 28907 voutt = fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v4( \ 28908 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt, \ 28909 di1=vshape[2], di2=vshape[1], di3=vshape[0], ds1=ssh[3], \ 28910 ds2=ssh[2], ds3=ssh[1], ds4=ssh[0], maxngridsin=sin.shape[1]) 28911 vout = voutt.transpose() 28912 elif len(vshape) == 3 and len(sN.shape) == 3: 28913 voutt = fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v3( \ 28914 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt, \ 28915 di1=vshape[2], di2=vshape[1], di3=vshape[0], ds1=ssh[2], \ 28916 ds2=ssh[1], ds3=ssh[0], maxngridsin=sin.shape[1]) 28917 vout = voutt.transpose() 28918 else: 28919 print errormsg 28920 print ' ' + fname + ": variable rank:", len(vshape), 'combined with', \ 28921 ' slice rank:', len(sN.shape), 'not ready!!' 28922 print ' available ones _______' 28923 print ' var_rank: 3 slice_rank: 3' 28924 print ' var_rank: 3 slice_rank: 4' 28925 quit(-1) 28926 28927 #print "Lluis Let's check !!", vout.shape 28928 #print 'i j k it:', i, j, k ,it, 'N:', sN[i,j,k], '<>', slicesvals 28929 #for iv in range(sN[i,j,k]): 28930 # print iv, 'c:', sin[:,iv,i,j,k], 'p:', sgpa[iv,i,j,k], 'v:', \ 28931 # varv[it,sin[0,iv,i,j,k],sin[1,iv,i,j,k]] 28932 #print 'min max mean mean2 std median count _______' 28933 #print vout[:,it,i,j,k] 28934 28935 ist = 0 28936 for statn in statns: 28937 newvar = onewvars[statn] 28938 rightvals = np.where(vout[6] == 0, gen.fillValueF, vout[ist]) 28939 newvar[:] = rightvals 28940 onewnc.sync() 28941 ist = ist + 1 28942 28943 onewnc.sync() 28924 28925 # Checking for size... 28926 Nvdims = len(vshape) 28927 varslc = {} 28928 if Nvdims == 3: 28929 vSIZE = np.prod(varv.shape[:])*8. 28930 if vSIZE > gen.oneGB*3.: 28931 print ' ' + infmsg 28932 print ' variable size: ', vSIZE / gen.oneGB,' larger than 3GB !!' 28933 print ' splitting time evolution to compute' 28934 vHSIZE = np.prod(varv.shape[1:3])*8. 28935 print ' var. horizontal size:', vHSIZE, 'dimt:', varv.shape[0] 28936 tsteps = gen.oneGB*3. / vHSIZE 28937 tslices = range(0,varv.shape[2],tsteps) 28938 print ' time-slices:', tklices 28939 for itt in tslcices[1,len(tslices)]: 28940 iit = tslices[itt-1] 28941 eet = tslices[itt] 28942 varslc[itt-1] = [slice(iit,eet), slice(0,vshape[1]), \ 28943 slice(0,vbshape[2])] 28944 else: 28945 tslices = range(0,varv.shape[0]) 28946 varslc[0] = [slice(0,vshape[0]), slice(0,vshape[1]), \ 28947 slice(0,vshape[2])] 28948 28949 for itt in varslc.keys(): 28950 slcv = varslc[itt] 28951 varv = ovar[tuple(slcv)] 28952 d0 = varv.shape[0] 28953 d1 = varv.shape[1] 28954 d2 = varv.shape[2] 28955 varvt = varv.transpose() 28956 if Nvdims == 3 and len(sN.shape) == 4: 28957 voutt=fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v4( \ 28958 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt, \ 28959 di1=d2, di2=d1, di3=d0, ds1=ssh[3], ds2=ssh[2], ds3=ssh[1], \ 28960 ds4=ssh[0], maxngridsin=sin.shape[1]) 28961 elif Nvdims == 3 and len(sN.shape) == 3: 28962 voutt=fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v3( \ 28963 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt, \ 28964 di1=d2, di2=d1, di3=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0], \ 28965 maxngridsin=sin.shape[1]) 28966 else: 28967 print errormsg 28968 print ' ' + fname + ": variable rank:", Nvdims, 'combined with', \ 28969 ' slice rank:', len(sN.shape), 'not ready!!' 28970 print ' available ones _______' 28971 print ' var_rank: 3 slice_rank: 3' 28972 print ' var_rank: 3 slice_rank: 4' 28973 quit(-1) 28974 28975 vout = voutt.transpose() 28976 ist = 0 28977 for statn in statns: 28978 newvar = onewvars[statn] 28979 rightvals = np.where(vout[6] == 0, gen.fillValueF, vout[ist]) 28980 newvar[tuple(slcv)] = rightvals 28981 onewnc.sync() 28982 ist = ist + 1 28983 onewnc.sync() 28984 #print "Lluis Let's check !!", vout.shape 28985 #print 'i j k it:', i, j, k ,it, 'N:', sN[i,j,k], '<>', slicesvals 28986 #for iv in range(sN[i,j,k]): 28987 # print iv, 'c:', sin[:,iv,i,j,k], 'p:', sgpa[iv,i,j,k], 'v:', \ 28988 # varv[it,sin[0,iv,i,j,k],sin[1,iv,i,j,k]] 28989 #print 'min max mean mean2 std median count _______' 28990 #print vout[:,it,i,j,k] 28944 28991 28945 28992 # Add global attributes
Note: See TracChangeset
for help on using the changeset viewer.
