[MITgcm-support] Sample pTracer using flt package
Martin Losch
Martin.Losch at awi.de
Thu Jun 21 04:55:37 EDT 2018
I would use the diagnostics package to save the tracer. Define a new diagnostic analogous to the TRACER01 …, (e.g. MYTRAC01) in ptracers_diagnostics_init and then add an appropriate CALL DIAGNOSTICS_FILL.
in the namelist you can control when in the timestep you want output.
timeaverages with frequency > 0
snapshots with frequency < 0
Note, that by default the time average is save at the end of the averaging interval (because it’s not possible to save it early, should be obvious), but the snapshot is save in the middle of the output interval. So if you have freqency = -86400, you’ll have output at 43200s + n*86400s. You can change it by setting timePhase (default frequencey/2) to your desired phase.
Martin
> On 21 Jun 2018, at 09:55, 钱钰坤 <qianyk at mail3.sysu.edu.cn> wrote:
>
> Hi Martin,
>
> I understand it much better now that we only pass a pointer to the S/R FLT_BILINEAR.
>
> I move all gradient-calculation codes from ptracer_integrate.F to flt_traj.F and add local
> fields of _RS xA/yA of one follow your suggestions. Now it works perfectly.
>
> Now I have online- and offline-sampled pTracer values and differences seem unacceptable
> for my purpose. After double-check the c-grid coordinates I used in offline way (may offset
> half grid length), I doubt that the instantaneous tracer outputs by diagnostics pkg (used for
> offline sampling) are somewhat different from the tracer state sampled online. One possible
> reason is that the model use AB integrator and has three timesteps (not sure if it is relevant).
> Another reason is that the output frequency is daily and the states may be defined exactly at
> the beginning of the day or the mid-day.
>
> In order to verify the differences, I just want to output the pTracer field just before Lagrangian
> sampling (calling flt_bilinear S/R) in flt_traj.F. I added the following codes:
> _BARRIER
> iRec = 1 + NINT( (myTime-startTime)/ PTRACERS_dumpFreq )
> CALL WRITE_REC_XYZ_RL('myPTracer',
> & pTracer(1-OLx,1-OLy,1,1,1,9),iRec,myIter,myThid)
> _BARRIER
> and make use of the predefined S/R write_rec_xyz_rl. First, this reported an error of 'myPTracer.data'
> file cannot be found. After I created an empty file of 'myPTracer.data', the code reported an error:
> forrtl: File too large
> forrtl: severe (38): error during write, unit 9, file myPTracer.data
> and the myPTracer.data is still 0 byte. I hope to get a single mdsio-like myPTracer.data binary file (not
> one file per tile). Is this error caused by not open a file or not flush the buffer to file? How to do it correctly?
> Are the _BARRIER lines necessary to synchronize all processes to generate a single file?
>
>
> ------------------
> Best regards
>
> Yu-Kun Qian (钱钰坤)
> Center for Monsoon and Environment Research
> Department of Atmospheric Sciences
> School of Environmental Science and Engineering
> Sun Yat-sen University
> No. 135 Xingang West Road, Haizhu District
> Guangzhou, 510275, P.R. China
> Tel; 020-84115227
> Email: qianyk at mail3.sysu.edu.cn
>
>
>
>
> ------------------ 原始邮件 ------------------
> 发件人: "Martin Losch"<Martin.Losch at awi.de>;
> 发送时间: 2018年6月14日(星期四) 晚上9:39
> 收件人: "MITgcm Support"<mitgcm-support at mitgcm.org>;
> 主题: Re: [MITgcm-support] Sample pTracer using flt package
>
> Hi Qian,
>
> I think I said before that I never used the flt-package, so I have no experience, how it actually works, but I’ll give ti a shot to answer you questions, see below:
> > On 14. Jun 2018, at 15:16, 钱钰坤 <qianyk at mail3.sysu.edu.cn> wrote:
> >
> > Hi Martin,
> >
> > Thanks for your help. I've compiled and run these codes successfully. However, I get all zeros after the run.
> >
> > 1) You have suggested to modify the Lagrangian sampling code as (see previous email below):
> > CALL FLT_BILINEAR (ix,jy,ss, pTracer(1-Olx,1-Oly,1,1,1, iTracer), kp,0,bi,bj,myThid)
> > I just wonder why this line is not like:
> > CALL FLT_BILINEAR (ix,jy,ss, pTracer(1-Olx,1-Oly,1,bi,bj, iTracer), kp,0,bi,bj,myThid)
> > The index of 1 rather than bi/bj seems to works fine for sampling the tracer field.
> in s/r flt_bilinear, the input “var” is defined as
> _RL var(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> so as a full 5D field, but you are passing ptracers which is 6D. In Fortran 77 you actually only pass the address of the first element of a field, and this is what I suggested by pTracer(1-Olx,1-Oly,1,1,1, iTracer), so s/r flt_bilinear “sees” the address of where to start and assumes a 5D field (the last index varies most slowly). I assume that your bi/bj = 1 anyway, so the other form should have worked (as long as you not using the multi-thread-code where bi/bj can be > 1), but better do it right the first time. Passing the indices bi,bj tell s/r flt_bilinear, where to operate on in the input file “var”. E.g. in pkg/gchem/gchem_calc_tendency.F, cfc11_forcing assumes 3D input fields and you have to pass the proper pTracer(1-Olx,1-Oly,1,bi,bj, iTracer) fields. In Fortran90 things would be more transparent.
> >
> > 2) the predefined subroutines GAD_GRAD_X/Y in GAD package return tracer gradient much
> > larger than expected. I looked into the code and found that these subroutines calculate the
> > product of gradient and face area of grid cell. So I may need to divide the face area to get
> > the right values of gradient. Is this correct? If so, I cannot just modify the subroutines because
> > they may be called in other places.
> Correct, it calculates a flux through the cell face, with area xA/yA
> >
> > 3) Subroutines GAD_GRAD_X/Y return non-zero values at the end of ptracer_integrate.F (this is
> > verified by outputting pTrGrdx/y after calling GAD_GRAD_X/Y subroutines), however,
> > in flt_traj.F, pTrGrdx/y becomes all zeros (this is verified by outputting pTrGrdx/y before calling
> > FLT_BILINEAR). Therefore, Lagrangian sampling gives all zero results. Does this related to any
> > execution sequence? Or any other possible reasons I cannot tell.
>
> how do you define your xA/yA in flt_traj? in ptracers_integrate they are computed in s/r calc_adv_flow. In your case, because you don’t want your gradients multiplied with the area you could just define local fields
> _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
> and set them to one. That would also solve your question (2).
>
> Martin
>
>
> >
> > Many thanks.
> >
> > ------------------
> > Best regards
> >
> > Yu-Kun Qian (钱钰坤)
> > Center for Monsoon and Environment Research
> > Department of Atmospheric Sciences
> > School of Environmental Science and Engineering
> > Sun Yat-sen University
> > No. 135 Xingang West Road, Haizhu District
> > Guangzhou, 510275, P.R. China
> > Tel; 020-84115227
> > Email: qianyk at mail3.sysu.edu.cn
> >
> >
> >
> >
> > ------------------ 原始邮件 ------------------
> > 发件人: "Martin Losch"<Martin.Losch at awi.de>;
> > 发送时间: 2018年6月14日(星期四) 晚上7:38
> > 收件人: "MITgcm Support"<mitgcm-support at mitgcm.org>;
> > 主题: Re: [MITgcm-support] Sample pTracer using flt package
> >
> > Looks good to me.
> >
> > M.
> >
> > > On 14. Jun 2018, at 05:23, 钱钰坤 <qianyk at mail3.sysu.edu.cn> wrote:
> > >
> > > Hi all,
> > >
> > > Still along this line, I want to sample pTracer along Lagrangian trajectory using flt
> > > package. Now I'm able to sample pTracer itself and I want to sample its gradient
> > > too.
> > >
> > > I add two arrays to store the gradient (vector) of pTracer in PTRACER_FIELDS.h:
> > > _RL pTrGrdx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> > > _RL pTrGrdy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
> > >
> > > Then I add the following lines at the end of ptracer_integrate.F to calculate the
> > > gradient of the ninth tracer, which makes use of the predefined subroutines of
> > > GAD_GRAD_X and GAD_GRAD_Y in GAD package:
> > > IF (iTracer .EQ. 9) THEN
> > > CALL GAD_GRAD_X(
> > > I bi,bj,k,xA,pTracer(1-Olx,1-Oly,1,bi,bj,9),
> > > O pTrGrdx(1-Olx,1-Oly,1,bi,bj),
> > > I myThid)
> > > CALL GAD_GRAD_Y(
> > > I bi,bj,k,yA,pTracer(1-Olx,1-Oly,1,bi,bj,9),
> > > O pTrGrdy(1-Olx,1-Oly,1,bi,bj),
> > > I myThid)
> > > ENDIF
> > >
> > > Finally, I modify the Lagrangian sampling code in flt_traj.F as:
> > > CALL FLT_BILINEAR (ix,jy,uu,pTrGrdx(1-Olx,1-Oly,1,1,1),
> > > I kp,1,bi,bj,myThid)
> > > CALL FLT_BILINEAR (ix,jy,vv,pTrGrdy(1-Olx,1-Oly,1,1,1),
> > > I kp,2,bi,bj,myThid)
> > > so that the gradients are stored in the u/v samples of particles.
> > >
> > > I just want to make sure that this is OK. I am aware of that gradX is defined
> > > at u-point and gradY is defined at v-point. So I guess this should be the best
> > > way to get online sampled tracer gradient.
> > >
> > >
> > > ------------------
> > > Best regards
> > >
> > > Yu-Kun Qian (钱钰坤)
> > > Center for Monsoon and Environment Research
> > > Department of Atmospheric Sciences
> > > School of Environmental Science and Engineering
> > > Sun Yat-sen University
> > > No. 135 Xingang West Road, Haizhu District
> > > Guangzhou, 510275, P.R. China
> > > Tel; 020-84115227
> > > Email: qianyk at mail3.sysu.edu.cn
> > >
> > >
> > >
> > >
> > > ------------------ 原始邮件 ------------------
> > > 发件人: "qianyk"<qianyk at mail3.sysu.edu.cn>;
> > > 发送时间: 2018年6月2日(星期六) 晚上7:12
> > > 收件人: "mitgcm-support"<mitgcm-support at mitgcm.org>;
> > > 主题: 回复： [MITgcm-support] Sample pTracer using flt package
> > >
> > > Hi Martin,
> > >
> > > Thanks very much for your suggestions. I can sample the value of ptracer online now.
> > >
> > > However, there are considerable differences between online-sampled and offline-sampled
> > > values. I've doubted about the results of my offline-sampling program and that is why I
> > > insist to get the online results to verify that. Since I have to also sample the gradient of
> > > the tracer, there are more questions follow this line.
> > >
> > > For a quick solution, I hope to compute the tracer gradient online and then sample the
> > > gradient field exactly the same way as I sample the tracer field (just passing the pointer
> > > of the gradient to the subroutine FLT_BILINEAR). However, I've no idea how and where
> > > to do this. I've found many diagnostics in available_diagnostics.log related to the tracer
> > > itself (UTRAC01, ADVrTr01, ...) which could be modified slightly to be the tracer gradient.
> > > I just wonder if this is the best way to do get Lagrangian sampling of tracer gradient.
> > >
> > > From a fundamental way, I just want to figure out the differences in the programs. The
> > > tracer field is output as MDSIO format. As MITgcm uses C-grid, should the position of the
> > > lower-left-corner value of the tracer be exactly zero, or be DXC (DYC), or be DXC/2 (DYC/2)?
> > > I think this tiny differences in position will cause considerable sampled tracer value if the
> > > gradient is large.
> > >
> > > I also set north and east wall boundaries so the domain is closed (the eastmost Nx and
> > > northmost Ny grid points are all set to the undefined value of 0). Could the Lagrangian
> > > particle get into the cell of these grid points or they just stay within [1 Nx-1]*[1 Ny-1]
> > > grid points? More specifically, what are the upper and lower bounds of position for a
> > > Lagrangian particle?
> > >
> > > Thanks again.
> > >
> > >
> > >
> > > ------------------
> > > Best regards
> > >
> > > Yu-Kun Qian (钱钰坤)
> > > Center for Monsoon and Environment Research
> > > Department of Atmospheric Sciences
> > > School of Environmental Science and Engineering
> > > Sun Yat-sen University
> > > No. 135 Xingang West Road, Haizhu District
> > > Guangzhou, 510275, P.R. China
> > > Tel; 020-84115227
> > > Email: qianyk at mail3.sysu.edu.cn
> > >
> > >
> > >
> > >
> > > ------------------ 原始邮件 ------------------
> > > 发件人: "Martin Losch"<Martin.Losch at awi.de>;
> > > 发送时间: 2018年6月1日(星期五) 晚上7:25
> > > 收件人: "MITgcm Support"<mitgcm-support at mitgcm.org>;
> > > 主题: Re: [MITgcm-support] Sample pTracer using flt package
> > >
> > > Hi Qian,
> > >
> > > I think you can do it that way.
> > > You need to include the ptracer variables at the beginning of the routine like this:
> > > #ifdef ALLOW_PTRACERS
> > > #include "PTRACERS_SIZE.h"
> > > #include “PTRACERS_FIELDS.h"
> > > #endif
> > >
> > > and it’s always a good idea to bracket pkg-specific code by CPP-flags:
> > > #ifdef ALLOW_PTRACERS
> > > … your code ...
> > > #endif
> > >
> > > I think you need to say
> > > CALL FLT_BILINEAR (ix,jy,ss, pTracer(1-Olx,1-Oly,1,1,1, iTracer), kp,0,bi,bj,myThid)
> > > to pass the approriate pointer to FLT_BILINEAR
> > >
> > > Martin
> > > > On 1. Jun 2018, at 12:19, 钱钰坤 <qianyk at mail3.sysu.edu.cn> wrote:
> > > >
> > > > Hi all,
> > > >
> > > > I've inserted an passive tracer (using ptracer package) and deployed
> > > > Lagrangian particles (using flt package) simultaneously into a dynamic flow.
> > > > I wanted to sample the tracer value along each Lagrangian float but I only
> > > > found the following five lines of code in flt_traj.F:
> > > >
> > > > CALL FLT_BILINEAR (ix,jy,uu,uVel, kp,1,bi,bj,myThid)
> > > > CALL FLT_BILINEAR (ix,jy,vv,vVel, kp,2,bi,bj,myThid)
> > > > CALL FLT_BILINEAR2D(ix,jy,pp,etaN, 0,bi,bj,myThid)
> > > > CALL FLT_BILINEAR (ix,jy,tt,theta, kp,0,bi,bj,myThid)
> > > > CALL FLT_BILINEAR (ix,jy,ss,salt, kp,0,bi,bj,myThid)
> > > >
> > > > which means the flt can only sample the basic flow states (u, v, eta, theta, and salt).
> > > >
> > > > I wonder if there is an easy way to achieve this by modifying the code slightly as:
> > > >
> > > > CALL FLT_BILINEAR (ix,jy,ss, pTracer(, , , , , iTracer), kp,0,bi,bj,myThid)
> > > >
> > > > and using the result to replace the salt which I do not care about?
> > > >
> > > > If yes, how to modify the #include or variable declarations at the beginning of flt_traj.F?
> > > >
> > > > Many thanks.
> > > >
> > > >
> > > > ------------------
> > > > Best regards
> > > >
> > > > Yu-Kun Qian (钱钰坤)
> > > > Center for Monsoon and Environment Research
> > > > Department of Atmospheric Sciences
> > > > School of Environmental Science and Engineering
> > > > Sun Yat-sen University
> > > > No. 135 Xingang West Road, Haizhu District
> > > > Guangzhou, 510275, P.R. China
> > > > Tel; 020-84115227
> > > > Email: qianyk at mail3.sysu.edu.cn
> > > >
> > > >
> > > > _______________________________________________
> > > > MITgcm-support mailing list
> > > > MITgcm-support at mitgcm.org
> > > > http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
> > >
> > > _______________________________________________
> > > MITgcm-support mailing list
> > > MITgcm-support at mitgcm.org
> > > http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
> > > _______________________________________________
> > > MITgcm-support mailing list
> > > MITgcm-support at mitgcm.org
> > > http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
> >
> > _______________________________________________
> > MITgcm-support mailing list
> > MITgcm-support at mitgcm.org
> > http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
> > _______________________________________________
> > MITgcm-support mailing list
> > MITgcm-support at mitgcm.org
> > http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mailman.mitgcm.org/mailman/listinfo/mitgcm-support
More information about the MITgcm-support
mailing list