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

lbc_lnk

lbc_lnk is done for all diagnostic variables in LOBSTER, e.g. trcbio.F90:

            CALL lbc_lnk( zw2d(:,:,jl),'T', 1. )

Q: do I need to do this in MY_TRC?

"tra" seems to represent the tracer concentration in the following time step, therefore I need to update tra instead of trn! The words "update the trend" are confusing to me (it seems the time rate of change is represented by tra. rather it should say update tra with the trend (or flux)).

tra, trn, trb for the next time step will be updated by TRP/trcnxt.F90, as well as the OBC and LBC are computed.

Actually, "tra" does represent the time rate of change in tracer concentration. For example PISCES/p4zflx.F90:

            ! compute the trend
            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) / fse3t(ji,jj,1)

where zfld - zflu represent the flux [mol/m2/s] which is divided by fse3t (vertical scale factor; thickness of first model layer??) [m]. Thus the division gives [mol/m3/s].

NaN in tracer.stat

I finally found the reason for segmentation fault (I think). It is because at least one of the newly added tracers (icechl, iceno3, and/or icenh4 in this case) has value of NaN, as seen in tracer.stat (which calculates sum of all 3d tracers, see ztrai in trcstp.F90):

         1  0.7883947556E+03
         2               NaN

This happens at the second time step. This means that either or both DMSPd or/and DMS becomes NaN in the second time step...

Actually this was not the cause of segmentation fault, as the run passed the second time step when key_my_trc_iceeco was excluded:

         1  0.7883947556E+03
         2               NaN
         3               NaN
         4               NaN
         5               NaN

Anyway, it is good to know that my 3d tracers are erraneous as well :(

time step for tracers in MY_TRC

PISCES/trcini_pisces.F90

rfact = rdttrc(1)

HH this should equal rdt (time step for dynamics) unless different tracer time step is given.

rfact2 = rfact / float (nrdttrc)

HH this would still be rfact unless nrdttrc (time step for biology) is different from tracers.

xstep = rfact2 / rday

HH rday = 86400 sec/day (defined in OPA_SRC/DOM/phycst.F90), so xstep gives the time step in days.. PISCES uses xstep when the rate constant is given in units of d-1. I will avoid using xstep and instead convert the rate constant to s-1 during the initialization as in FABM.