2016年10月6日木曜日

generating meteo file for amundsen


for month in `seq 1 12`;do wget --content-disposition "http://climate.weather.gc.ca/climate_data/bulk_data_e.html?format=csv&stationID=1633&Year=2008&Month=${month}&Day=14&timeframe=1&submit= Download+Data" ;done
Other stations in the Amundsen Gulf area:

2016年10月5日水曜日

When all three my_trc components are added, the tracer.stat increases too rapidly:

tail tracer.stat
         1  0.7895863787E+03
         2  0.3124101036E+04
         3  0.1067765821E+13
         4  0.2180091013E+27
         5  0.8173235638E+41

This still occurs even I remove the icedms component:

         1  0.7895863787E+03
         2  0.3056491046E+04
         3  0.1034367164E+13
         4  0.2180090884E+27
         5  0.8173235638E+41

When I removed the chunk of code about ice ablation/accretion, the tracer.stat became normal:

         1  0.7895862638E+03
         2  0.7895870356E+03
         3  0.7913792876E+03
         4  0.7936409329E+03
         5  0.7949360900E+03

Thus, the problem is here:

         ! Lateral accretion/ablation
        IF(frld(ji,jj)-pfrld(ji,jj) > 0.) THEN ! In case of accretion, newly formed ice will have a tracer concentration of an underlying seawater value.
          icechl(ji,jj) = (1.-pfrld(ji,jj))/(1.-frld(ji,jj)+epsilon10)*icechl(ji,jj) + ((1.-frld(ji,jj))-(1.-pfrld(ji,jj)))/(1.-frld(ji,jj)+epsilon10)*trn(ji,jj,1,jpdch)
          iceno3(ji,jj) = (1.-pfrld(ji,jj))/(1.-frld(ji,jj)+epsilon10)*iceno3(ji,jj) + ((1.-frld(ji,jj))-(1.-pfrld(ji,jj)))/(1.-frld(ji,jj)+epsilon10)*trn(ji,jj,1,jpno3)
          icenh4(ji,jj) = (1.-pfrld(ji,jj))/(1.-frld(ji,jj)+epsilon10)*icenh4(ji,jj) + ((1.-frld(ji,jj))-(1.-pfrld(ji,jj)))/(1.-frld(ji,jj)+epsilon10)*trn(ji,jj,1,jpnh4)
        END IF

2016年10月4日火曜日

interesting segfault

I had a segmentation fault which is due to NaN. However, this was interesting because although the IF statement was always false, it gave me segfault because NaN would result if the IF statement was true.. 

        write(numout,*) 'hhpass: nancheck ',(1.-pfrld(ji,jj))/(1.-frld(ji,jj)+epsilon30)
        IF(frld(ji,jj)-pfrld(ji,jj) > 0.) THEN ! In case of accretion, newly formed ice will have a tracer concentration of an underlying seawater value.
          write(numout,*) 'hhpass: inside if statement'
          icechl(ji,jj) = icechl(ji,jj)/(trn(ji,jj,1,jpdch)+epsilon30)!((1.-pfrld(ji,jj))/(1.-frld(ji,jj)+epsilon30)*icechl(ji,jj) + ((1.-frld(ji,jj))-(1.-pfrld(ji,jj)))/(1.-frld(ji,jj)+epsilon30)*trn(ji,jj,1,jpdch)
        END IF


 Update: This seems to have to do with floating points.

Writing out the output for the variables in dinomenator and numerator show a particular floating point (i.e. 14):

 hhpass: pfrld    1.00000000000000
 hhpass: frld    1.00000000000000





In the example above, I tried avoid the division by zero by introducing epsilon number 30 (i.e. 1.e-30), but the division produced NaN:

hhpass: nancheck    NaN

After this, I instead introduced a epsilon number that is smaller than the floating point of the two variables (i.e. 4). Then the division produced zero:

 hhpass: nancheck   0.000000000000000E+000


Epsilon30 worked for biogeochemical tracers. But not for ice physical tracers (probably a specic floating point number is defined in LIM_SRC_2 for these tracers?)

2016年10月3日月曜日

NaN in tracer.stat --> found the problem!

I used to see NaN in the second time step of tracer.stat. In the first time step you dont see Nan because it is the initial condition.

In the second time step, seawater NO3 was NaN because it was exchanging the NaN value coming from the sea ice NO3, which was due to the lack of epsilon number.

For example, the sink for sea ice NO3 due to ice algal uptake is calculated based on the ratio of NO3 and all other available nitrogen (i.e. NO3 and NH4). When NO3 and NH4 are initially set to 0 in the model (for now), the denominator becomes zero, therefore the sink = NaN, therefore the ice NO3 in the next time step as well as the oceanic NO3 became NaN.

So I fixed this by adding a very small epsilon number to avoid this division by zero.

interpolate seawifs onto NAA grid

1) Download the monthly climatology data:
wget http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19972442010273.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19972742010304.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19973052010334.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19973352010365.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19980012010031.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19980322010059.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19980602010090.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19980912010120.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19981212010151.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19981522010181.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19981822010212.L3m_MC_CHL_chlor_a_9km.nc http://oceandata.sci.gsfc.nasa.gov/cgi/getfile/S19982132010243.L3m_MC_CHL_chlor_a_9km.nc

2) Rename the data so that they will appear in the order of months (i.e. from Jan. to Dec.):
for i in $(ls *_MC_*.nc); do `echo "cp $i tmp${i:5}"`;done

ls tmp*.nc
tmp0012010031.L3m_MC_CHL_chlor_a_9km.nc
tmp0322010059.L3m_MC_CHL_chlor_a_9km.nc
tmp0602010090.L3m_MC_CHL_chlor_a_9km.nc
tmp0912010120.L3m_MC_CHL_chlor_a_9km.nc
tmp1212010151.L3m_MC_CHL_chlor_a_9km.nc
tmp1522010181.L3m_MC_CHL_chlor_a_9km.nc
tmp1822010212.L3m_MC_CHL_chlor_a_9km.nc
tmp2132010243.L3m_MC_CHL_chlor_a_9km.nc
tmp2442010273.L3m_MC_CHL_chlor_a_9km.nc
tmp2742010304.L3m_MC_CHL_chlor_a_9km.nc
tmp3052010334.L3m_MC_CHL_chlor_a_9km.nc
tmp3352010365.L3m_MC_CHL_chlor_a_9km.nc

3) Create a new variable called time to each data:
for i in $(ls tmp*9km.nc); do ncap2 -O -s 'time=array(1.,1.,1.)' $i vartime_$i; done

4) Create an unlimited (RECORD) dimension called "time" to each data:
 for i in $(ls vartime_*.nc); do ncecat -u time $i dimtime_$i; done

5) Concatenate the files:
ncrcat -h dimtime_*.nc seawifs_mon_clim.nc

6) Finally use SOSIE to interpolate onto NAA grid:
~/sosie-2.6.4/bin/sosie.x -f namelist.seawifs2naa