[MITgcm-support] passive tracer restoring
Martin Losch
Martin.Losch at awi.de
Fri Jul 28 03:27:46 EDT 2017
Hi Qian,
if the range valid tracer values is [0,someMax], then the function should return the correct values, i.e. zero for the minimum. I am not sure what you are trying to do. I can think of two options
(1) you don’t want to take into account tracer values of exactly 0 (which does not make too much sense to me), then your code does what you want.
(2) you don’t want the zero values on dry cells to contribute (they are all zero), but zero tracer values in the model interior (wet cells) are OK; in this case you should use maskC(I,J,k,bi,bj) instead of ptracers in the if-statement. That’s how the the monitor package does it. If you don’t want to take the values on the open boundaries into account either, then you need to use maskInC(i,j,bi,bj)*maskC(bi,bj,k,bi,bj) instead of maskC(I,J,k,bi,bj)
Martin
> On 27. Jul 2017, at 14:28, 钱钰坤 <qianyk at mail3.sysu.edu.cn> wrote:
>
> Hi Martin,
>
> Yes, I should put the call of the S/R outside the bi/bj loop.
>
> As for the ptracer = 0 problem, there are two points.
>
> 1) the theMax and theMin values are, in general, initialized by the possible minimum and maximum values respectively, in order to get the extrema. Now they are initialized by +/-1. _d 23, which may not be the maximum/minimum possible number for both real4 (order of magnitude ~38) or real8 (order of magnitude ~308) types. For a tracer like the streamfunction, the order of magnitude may exceed this, and then the extrema returned may be incorrect. Are there any predefined constants that can be used to initialize the max/min values without worrying about the real4 or real8 type?
>
> 2) the if statement for ptracer = 0 is added because without it I always get zero for theMin, while theMax is OK. Since there is no _GLOBAL_MIN_RL, you told me a solution of finding the max of minus theMin. In my case, ptracer are all positive. Once taking the minus, they all be negative. At the boundary there are zeros in ptracer field. So MAX S/R finally gets 0 and return it. If adding a IF statement is not the best, then can you suggest another way to fix it?
>
> 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月27日(星期四) 晚上6:28
> 收件人: "MITgcm Support"<mitgcm-support at mitgcm.org>;
> 主题: Re: [MITgcm-support] passive tracer restoring
>
> Hi Qian,
>
> looks good except that you use the bi/bj-loops twice. Because you want to do global max, you need to have call your new routine after the end of the bi/bj loops and do the bi/bj loops inside the routine, see attached code:
>
> C-- End of 2nd bi,bj loop
> ENDDO
> ENDDO
>
> C miniufo added: call regularization
> #ifdef ALLOW_PTRACERS
> IF (usePTRACERS)
> & CALL PTRACERS_REGULARIZE( myTime, myIter, myThid )
> #endif
>
> (1) I am not sure why you need to have this if statement for ptracers = 0 if you want to allow negative values at the same time, but that’s beyond my comprehension of this parameterization. If you jus want to mask land values it would make more sense to use maskC instead of ptracers.
> (2) In the attached code I added the k-loop (inside the bi/bj-loop) to compute max/min of 3D fields. As long as you run only with one level, you can use the code as it is (I did not check it even for compilation errors). [If you want min/max of 2D slices, you’ll have to move the k-loop outside of the bj-bj loops, basically just one big k-loop for the entire routine, but I guess that’s not what you want]
> (3) halos are for parallelization only. Check out the manual for an explanation (Software Architecture: <http://mitgcm.org/public/r2_manual/latest/online_documents/node181.html>). maskInC is 2D (only surface) and excludes points along open boundaries (“In” is "interior of the domain without open boundaries). If you want to use it on a 3D field, you need to use maskInC(i,j,bi,bj)*maskC(bi,bj,k,bi,bj)
>
> Martin
>
>
>
>
>
> > On 27. Jul 2017, at 10:32, 钱钰坤 <qianyk at mail3.sysu.edu.cn> wrote:
> >
> > Oops, missing attached files 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
> >
> >
> >
> >
> > ------------------ 原始邮件 ------------------
> > 发件人: "qianyk"<qianyk at mail3.sysu.edu.cn>;
> > 发送时间: 2017年7月27日(星期四) 下午4:26
> > 收件人: "mitgcm-support"<mitgcm-support at mitgcm.org>;
> > 主题: Re: Re: [MITgcm-support] passive tracer restoring
> >
> > 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
> > <ptracers_regularize.F><tracers_correction_step.F>_______________________________________________
> > 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
More information about the MITgcm-support
mailing list