[MITgcm-support] passive tracer restoring
钱钰坤
qianyk at mail3.sysu.edu.cn
Thu Jul 27 04:26:24 EDT 2017
Hi Martin,
Sorry for a bit delay as I spent some time to give a try. Thanks very much for your code snippet.
As you suggested, I added a S/R PTRACERS_REGULARIZATION.F (using PTRACERS_INTEGRATE.F as template) in the code directory
and then make a call at the end of TRACERS_CORRECTION_STEP.F. I've attached both files so that you may find some silly bugs or
redundant statements (e.g., #INCLUDE).
There are still some problems when I run these codes:
1) The theMin value is always zero. I think this is related to the default undefined value of zero used by the model. So I added an IF as:
theMax = -1. _d 23
theMin = 1. _d 23
DO bj = myByLo(myThid),myByHi(myThid)
DO bi = myBxLo(myThid),myBxHi(myThid)
DO j = 1, sNy
DO i = 1, SNx
IF (pTracer(i,j,k,bi,bj,iTracer) .NE. 0. _d 0) THEN
theMax = MAX(pTracer(i,j,k,bi,bj,iTracer), theMax)
theMin = MIN(pTracer(i,j,k,bi,bj,iTracer), theMin)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
so that the zero values of the ptracer do not enter the MAX/MIN funcitons. Besides, the theMax is initialized with negative theMin
(rather than 0 as you suggest), in case tracer values are all negative.
2) Currently, I running the model with only one level. In the future, I may use a multi-level setup. How to code the k-loop if I want to
find the max/min value of a 3D tracer field? Is it OK just adding a k-loop outside the bi/bj loop?
3) I don't get the concept of halo regions. Are they similar to isolated lakes from the ocean? As you mentioned you only use interior
of the domain, do you mean that the maskInC is preferred (rather than maskC)?
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>;
发送时间: 2017年7月20日(星期四) 下午4:17
收件人: "MITgcm Support"<mitgcm-support at mitgcm.org>;
主题: Re: [MITgcm-support] passive tracer restoring
Hi,
so your method is not really a restoring method, but a regularization, isn’t it?
Currently, you do this only for the surface layer (according to your code). If this is what you want, I suggest that you write the min/max yourself and don’t use the expensive mon_calc_stats_rl, which will compute the min/max of the full 3D field. Also by placing it in to ptracers_integrate, you repeat the computation for every process and call the global reductions from within the bi/bi loops. This will only work of you use one tile per MPI process (which you probably do).
We have functions for the global_max reductions. I would also use only the interior of the domain (without the overlaps or halos, who knows what’s in there at this stage, at best copies of some interior values, update of the halos will be done afterward in the blocking_exchange routines) and I would put the code somewhere outside the bi/bj-loops, e.g. at the end of tracers_correction_step.F; there I suggest to put this code (to be really clean, put it into a subroutine PTRACERS_REGULARIZE or similar):
C need to include PTRACERS.h PTRACERS_SIZE.h and define the appropriate
C variables here
DO iTracer = 1, PTRACERS_numInUse
k = 1
theMax = 0. _d 0
theMin = 1. _d 23
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1,sNy
DO i=1,sNx
theMax = MAX(pTracer(i,j,k,bi,bj,iTracer),theMax)
theMin = MIN(pTracer(i,j,k,bi,bj,iTracer),theMin)
ENDDO
ENDDO
ENDDO
C There is no _GLOBAL_MIN_RL( theMin, myThid )
theMin = -theMin
_GLOBAL_MAX_RL( theMin, myThid )
theMin = -theMin
_GLOBAL_MAX_RL( theMax, myThid )
C now you have the max/min values for layer k
maxInit = ...
minInit = ...
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1,sNy
DO i=1,sNx
pTracer(i,j,k,bi,bj,iTracer) = ( minInit +
& ( maxInit - minInit ) * (
& pTracer(i,j,k,bi,bj,iTracer) - theMin
& ) / ( theMax - theMin ) ) * maskInC(i,j,bi,bj)
&
ENDDO
ENDDO
ENDDO
ENDDO
C tracer loop
ENDDO
I hope I am making sense,
Martin
> On 20. Jul 2017, at 04:34, 钱钰坤 <qianyk at mail3.sysu.edu.cn> wrote:
>
> Hi Martin,
>
> Thanks very much for your suggestion. I agree with you that it is much simpler not to involve the RBCS package.
>
> After reading the source code of the ptracer pkg, I've added the restoring code at the end of the ptracers_integrate.F (see the attached figure), maybe not the best.
>
> The main problem is that the restoring code requires the maximum and minimum values of the tracer field at each timestep. This is a reduction operation in parallel environment and is not simple for me to write a subroutine that returning the max and min. Fortunately, I found the existing subroutine MON_CALC_STATS_RL which is in monitor pkg and summarizes statistics for any given fields (actually, I only need max and min values while this subroutine provides more than this and thus could slow down the integration if it is called each timestep).
>
> Now the code works fine but I still not sure about:
> 1) the upper and lower bound of the index for the loop (i,j) is correct (1-OLy to sNy+OLy for j while 1-OLx to sNx+OLx for i);
> 2) the use of maskInC or maskC to maskout the topo (I find maskInC works fine).
>
> 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>;
> 发送时间: 2017年7月18日(星期二) 晚上7:56
> 收件人: "MITgcm Support"<mitgcm-support at mitgcm.org>;
> 主题: Re: [MITgcm-support] passive tracer restoring
>
> Hi,
>
> did you get an answer?
>
> You could modify the RBCS package to do what you want to do, but maybe simpler to edit ptracers_apply_forcing.F, where you can add any passive tracer tendency. You need to reformulate your method so that translates into a tendency (tracerunits/time) and add it to gPtracer in an extra loop. This routine is called for every vertical level (k), so a 2D loop sufficies:
>
> DO j=jMin,jMax
> DO i=iMin,iMax
> gPtracer(i,j) = gPtracer(i,j) + yourTendency(i,j)
> ENDDO
> ENDDO
>
> Martin
>
> > On 7. Jul 2017, at 05:25, 钱钰坤 <qianyk at mail3.sysu.edu.cn> wrote:
> >
> > Hi all,
> >
> > I've released a passive tracer q in a specific setup and let the tracer evolve under the dynamic flows without any sources.
> >
> > Due to the numerical diffusion induced by advection schemes (no explicit diffusion is set), the tracer variances decay slowly with time. Therefore, for a long-time run, the global averaged tracer variance (or overall gradient estimated roughly by qmax-qmin) will decrease all the time. To prevent this, I need the restoring method by Alen and Nakamura (2003, JAS):
> >
> > q*=-1+2(q-qmin)/(qmax-qmin)
> >
> > After this restoring (and normalization) at each time step, the tracer maximum and minimum values will remain between -1 and 1, and do not decay with time.
> >
> > The problem is where to add this restoring code (or which source file to add to) every time step, so that it is a simple and easy task.
> >
> > I've thought to use RBCS but RBCS seems to suitable for prescribed tracer field and the restoring formula here is tracer dependent.
> >
> > Thanks in advance.
> >
> > Reference:
> > Allen, D. R. and N. Nakamura, 2003: Tracer Equivalent Latitude: A Diagnostic Tool for Isentropic Transport Studies. J. Atmos. Sci., 60, 287-304.
> >
> >
> > ------------------
> > 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://mitgcm.org/mailman/listinfo/mitgcm-support
>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
_______________________________________________
MITgcm-support mailing list
MITgcm-support at mitgcm.org
http://mitgcm.org/mailman/listinfo/mitgcm-support
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20170727/6d38bd0b/attachment-0001.htm>
More information about the MITgcm-support
mailing list