[ncl-talk] rho_mwjf

[Chathurika]

Dear NCL users,




I am really appreciate if you can show me what I is wrong with my script. As you know the southern ocean density is around 26-27 as follows at each depth level.




(0) 26.34497

(1) 26.3656

(2) 26.40527

(3) 26.45532

(4) 26.5199

(5) 26.58557

(6) 26.65369

(7) 26.71924

(8) 26.77869

(9) 26.83704

(10) 26.89197

(11) 26.94214

(12) 26.99036

(13) 27.0343

(14) 27.0769

(15) 27.11511

(16) 27.15247

(17) 27.18738

(18) 27.22571

(19) 27.26978

(20) 27.32153

(21) 27.39014

(22) 27.46655

(23) 27.52747

(24) 27.57336

(25) 27.61084

(26) 27.65247

(27) 27.6958

(28) 27.73132

(29) 27.75598

(30) 27.77417

(31) 27.78784

(32) 27.79846

(33) 27.80762

(34) 27.81592

(35) 27.82202

(36) 27.82605

(37) 27.82739

(38) 27.82776

(39) 27.8324




However, when I calculate the density at each depth levels using rho_mwjf (script is attached below), It calculate density which is not reasonable, shows;




(0) 26.3907

(1) 26.4604

(2) 26.54432

(3) 26.6379

(4) 26.74557

(5) 26.85472

(6) 26.97386

(7) 27.10243

(8) 27.24155

(9) 27.40277

(10) 27.58386

(11) 27.78384

(12) 28.00518

(13) 28.24571

(14) 28.5084

(15) 28.79022

(16) 29.09488

(17) 29.43268

(18) 29.82106

(19) 30.26328

(20) 30.76077

(21) 31.32516

(22) 31.94666

(23) 32.59835

(24) 33.27908

(25) 33.99501

(26) 34.78358

(27) 35.64154

(28) 36.53573

(29) 37.46336

(30) 38.45152

(31) 39.54657

(32) 40.74771

(33) 42.0564

(34) 43.53943

(35) 45.23687

(36) 47.14139

(37) 49.24323

(38) 51.43494

(39) 53.81633




depth from 6 m to 5720 m, 40 levels





Could please anyone say what am I doing here wrong? if this is not the right way to calculate density using temperature and salinity what will be the other way? Many many thanks. Please be kind to reply me.




The script: 




diri=".../salinity/hist/"
sfile = "salinity.nc" 


f = addfile(diri+sfile,"r")
vname=getfilevarnames(f)
print(vname)


  salt           = f->so(:,:,{-25:-70},:)
  printVarSummary(salt)
  printMinMax(salt, False)


 s_range = dim_avg_n_Wrap(salt, 0)
 printVarSummary(s_range)


 delete(salt)


diri=".../temperature/hist/"
sfile = "temperature.nc" 


f = addfile(diri+sfile,"r")
vname=getfilevarnames(f)
print(vname)


  temp           = f->thetao(:,:,{-25:-70},:)
  printVarSummary(temp)
  printMinMax(temp, False)


 t_range = dim_avg_n_Wrap(temp, 0)
 printVarSummary(t_range)


 delete(temp)


depth = (/6, 17,27,37,47,57,68.5,82.5,100,122.5,150,182.5,220\
   ,262.5,310,362.5,420,485,560,645,740,845,960,1085,1220,1365\
   ,1525,1700,1885,2080,2290,2525,2785,3070,3395,3770,4195,4670,5170,5720/) 


ndim = dimsizes(t_range)


nlev = ndim(0)
nlat = ndim(1)
nlon = ndim(2)




rho = new((/nlev,nlat,nlon/), typeof(t_range))


do i = 0, nlev-1
    rho(i,:,:) = rho_mwjf(t_range(i,:,:),s_range(i,:,:),depth(i))
   
end do


copy_VarCoords(t_range, rho)
rho = 1000.*(rho-1.) 


printVarSummary(rho)
printMinMax(rho, True)


abc = dim_avg_n_Wrap(rho, (/1,2/))
printVarSummary(abc)


print(abc)





Thank you and best regards,

Chathu



Wickramage Chathurika Hemamali
Msc in Physical Oceanography

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220115/a10f0651/attachment.html>