<div dir="ltr">Hello all,<div><br></div><div>I have made a .ctl file from .nc file using CDO command which shows following description-</div><div><br></div><div><div>* Generated by CDO operator gradsdes</div><div>*</div><div>DSET  ^<a href="http://merge_rf79_05_act_0lag.nc">merge_rf79_05_act_0lag.nc</a></div><div>DTYPE  netCDF</div><div>XDEF 720 LINEAR 0.000000 0.500000</div><div>YDEF 361 LINEAR -90.000000 0.500000</div><div>ZDEF 7 LEVELS 2 3 5 7 8.5 9.25 10 </div><div>TDEF 164 LINEAR 12:00Z03aug1979 1dy</div><div>TITLE  <a href="http://merge_rf79_05_act_0lag.nc">merge_rf79_05_act_0lag.nc</a>  720x361 grid</div><div>OPTIONS yrev zrev</div><div>UNDEF  -32767</div><div>VARS  4</div><div>v                  7  t,z,y,x  V component of wind  [m s**-1]</div><div>u                  7  t,z,y,x  U component of wind  [m s**-1]</div><div>t                  7  t,z,y,x  Temperature  [K]</div><div>w                  7  t,z,y,x  Vertical velocity  [Pa s**-1]</div><div>ENDVARS</div></div><div><br></div><div>Using this .ctl file, i am trying to write a .ctl file containing variables for my diabatic heating calculation but the new .ctl file is not opening and it gives following error-</div><div><br></div><div><div>Scanning description file:  thermo.ctl</div><div>Open Error:  Looking for &quot;endvars&quot;, found &quot;res   7  99  residuo of diabatic heating  (q)  [k/dia] &quot; instead.</div><div>  --&gt; The invalid description file record is:</div><div>  --&gt; res   7  99  residuo of diabatic heating  (Q)  [K/dia]</div><div>  The data file was not opened.</div></div><div>I am using the following script to calculate the diabatic heating variable in .ctl file-</div><div><div>&#39;reinit&#39;</div><div>&#39;open merge_rf79_05_act_0lag.ctl&#39;</div><div><br></div><div>***********************************************************</div><div>&#39;q file&#39;</div><div>rc=sublin(result,5)</div><div>nt=subwrd(rc,12)</div><div>nl=subwrd(rc,9)</div><div>nx=subwrd(rc,3)</div><div>ny=subwrd(rc,6)</div><div>rc=sublin(result,6)</div><div>&#39;q time&#39;</div><div>td=subwrd(result,3)</div><div>ano=substr (td, 9, 4)</div><div><br></div><div>***********************************************************</div><div><br></div><div>&#39;set fwrite thermo.bin&#39;</div><div>&#39;set gxout fwrite&#39;</div><div>&#39;set x 1 720 &#39;</div><div>&#39;set y 1 361 &#39;</div><div><br></div><div>* Note: Repeat 1000 and 100 to permit program diferentiate z levels.</div><div><br></div><div>string=&#39; 10 10 9.25 8.5 7 5 3 2 2 &#39;</div><div><br></div><div>*aa      --&gt; m</div><div>*Rgas    --&gt; 287.05 J/Kg K</div><div>*Cp      --&gt; 1004   J/Kg K</div><div>*dt      --&gt; s</div><div>*T       --&gt; K</div><div>*u       --&gt; m/s</div><div>*v       --&gt; m/s</div><div>*P       --&gt; Pa</div><div>*W       --&gt; Pa/s</div><div><br></div><div>&#39;define aa=6.37e6&#39;</div><div>&#39;define pi=2*asin(1)&#39;</div><div>&#39;define rd=pi/180&#39;</div><div>&#39;define clat=cos(lat*rd)&#39;</div><div>&#39;define dx=cdiff(lon,x)*rd&#39;</div><div>&#39;define dy=cdiff(lat,y)*rd&#39;</div><div>&#39;define kp=287.05/1004&#39;</div><div>&#39;define dt=3600*6*2&#39;</div><div><br></div><div>&#39;k = 1&#39;</div><div>&#39;l = 2&#39;</div><div><br></div><div>while (k &lt;= 164)</div><div><br></div><div>while (l &lt;= 7+1)</div><div><br></div><div>&#39;set t &#39;k</div><div><br></div><div>l1=l-1</div><div>l2=l+1</div><div><br></div><div>nivel  = subwrd(string,l)</div><div>nivel1 = subwrd(string,11)</div><div>nivel2 = subwrd(string,l2)</div><div><br></div><div>* As vertical velocity is in units of Pa/sec and the levels (in CTL and nc data) </div><div>* are in hPa, we must multiply the levels (hPa) by 100 to convert it in Pa units </div><div>* (the same unit of vertical velocity)</div><div><br></div><div>&#39;define p1=&#39;nivel1&#39;*100&#39;</div><div>&#39;define p2=&#39;nivel2&#39;*100&#39;</div><div><br></div><div>say nivel&#39; &#39;nivel1&#39; &#39;nivel2&#39; &#39;k&#39; &#39;l</div><div><br></div><div>&#39;set lev &#39;nivel</div><div><br></div><div>* Computing Temperature tendency term</div><div>*_____________ dT/dt _______________*</div><div>* Units must be in K/s</div><div><br></div><div>k1=k-1</div><div>k2=k+1</div><div><br></div><div><br></div><div>if(k!=1 | k!=nt)</div><div>&#39;define ttend = (T(t=&#39;k2&#39;)-T(t=&#39;k1&#39;))/dt&#39;</div><div>endif</div><div><br></div><div>if(k=1)</div><div>&#39;define ttend = (T(t=&#39;k2&#39;)-T(t=&#39;k&#39;))/(dt/2)&#39;</div><div>say &#39;k=1&#39;</div><div>endif</div><div><br></div><div><br></div><div>if(k=nt)</div><div>&#39;define ttend = (T(t=&#39;k&#39;)-T(t=&#39;k1&#39;))/(dt/2)&#39;</div><div>say &#39;k=nt&#39;</div><div>endif</div><div><br></div><div>* Computing Temperature horizontal advection</div><div>* in hespherical coordinates</div><div>*_____________ V Grad(T) _______________*</div><div>* Units must be in K/s</div><div><br></div><div>&#39;define advh = ((u*cdiff(T,x))/(clat*dx) + (v*cdiff(T,y))/dy)/aa&#39;</div><div><br></div><div>* Computing Temperature vertical advection</div><div>* in hespherical coordinates</div><div>*_________ w Grad(T) or wdT/dp __________*</div><div>* Units must be in K/s</div><div><br></div><div>&#39;define advv = w*(t(lev=&#39;nivel2&#39;)-t(lev=&#39;nivel1&#39;))/(p2-p1)&#39;</div><div><br></div><div>* Computing alfa term</div><div>* in hespherical coordinates</div><div>* _____RTW/CpP _____________*</div><div>* Units must be in K/s</div><div><br></div><div>&#39;define qd   = -kp *  (w*(t))/(lev*100)&#39;</div><div><br></div><div>* Computing ther esidual term or the diabatic term</div><div>* in hespherical coordinates</div><div>* _____ Q _____________*</div><div>* Units must be in K/s</div><div><br></div><div>&#39;define res   = ttend + advh + advv + qd&#39;</div><div><br></div><div>* Recording the diabatic term in units of K/day</div><div><br></div><div>&#39;d ttend*86400&#39;</div><div>&#39;d advh*86400&#39;</div><div>&#39;d advv*86400&#39;</div><div>&#39;d qd*86400&#39;</div><div>&#39;d res*86400&#39;</div><div><br></div><div><br></div><div><br></div><div>k = k</div><div>l = l + 1</div><div><br></div><div>endwhile</div><div><br></div><div>l = 2</div><div>k = k + 1</div><div><br></div><div>endwhile</div><div><br></div><div>&#39;disable fwrite&#39;</div><div><br></div><div>**********************************************************</div><div>write (thermo.ctl,&#39;dset thermo.bin&#39;)</div><div>write (thermo.ctl,&#39;undef  -32767&#39;)</div><div>write (thermo.ctl,&#39;title <a href="http://merge_rf79_05_act_0lag.nc">merge_rf79_05_act_0lag.nc</a>  720x361 grid&#39;)</div><div>write (thermo.ctl,&quot;xdef &quot;720&quot; linear    0.0 &quot;0.5&quot;&quot;)</div><div>write (thermo.ctl,&quot;ydef &quot;361&quot; linear  -90.0 &quot;0.5&quot;&quot;)</div><div>write (thermo.ctl,&quot;zdef &quot;7&quot; levels  &quot;2&quot; &quot;3&quot; &quot;5&quot; &quot;7&quot; &quot;8.5&quot; &quot;9.25&quot; &quot;10&quot; &quot;)</div><div>write (thermo.ctl,&quot;tdef &quot;164&quot; linear &quot;00Z03aug1979&quot; &quot;12hr&quot; &quot;   )</div><div>write (thermo.ctl,&#39;vars  4&#39;)</div><div>write (thermo.ctl,&#39;ttend 7  99  temperature tendency  [K/dia] &#39;)</div><div>write (thermo.ctl,&#39;advh  7  99  horiz. temperature advection  [K/dia] &#39;)</div><div>write (thermo.ctl,&#39;advv  7  99  vertical temperature advection  [K/dia] &#39;)</div><div>write (thermo.ctl,&#39;qd    7  99                                  [K/dia] &#39;)</div><div>write (thermo.ctl,&quot;res   &quot;7&quot;  99  residuo of diabatic heating  (Q)  [K/dia] &quot;)</div><div>write (thermo.ctl,&#39;endvars&#39;)</div><div>**********************************************************</div></div><div><br></div><div>Sorry for the long note. Any help will be much appreciated.</div><div><br></div><div>Thanks</div><div>Praveen</div><div><br></div></div>