MPDATA and MASKING

[adil]

Dear Modelers

I performed some simulations with a modified GRAV_ADJ case with the latest source code which I downloaded a few weeks ago, the model resulted in really different solutions when my rectangular basin is closed with MASKING option instead of the WALL CPPFLAGS, actually when the masking is activated mpdata scheme finds really lower salinity values than the minumum of the initial field (without any forcing), this difference increases when the grid size is increased. When MPDATA is switched with TS_U3HADVECTION and TS_SVADVECTION schemes, masking and wall definitions don't differ much.

I did lots of tests to resolve this situation, finally I thought that the MPDATA advection scheme may have a little bug when used with the masking option. After some searching in the code, I concluded that the problem is; mpdata_adiff.F's supressing false oscillations by limiting the transport fluxes without skipping the masked points, it is better to consider the mask array while calculating the Max and Min of the advected tracer. My solution is (with line numbers) modifying the mpdata_adiff.F as:

adding a mask_min array

Code: Select all

143       real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: mask_min

setting Ta and Tadv to 1.0 at the masked points and setting mask_min array to 100.0 at land (or another high number which depends on your advected tracer) 1.0 at sea. I also set beta_up and beta_dn to 2.0 for the masked points.

Code: Select all

870 !  Scale advected tracer to Tunits.
 871 !
 872       DO k=1,N(ng)
 873         DO j=-1+J_RANGE+1
 874           DO i=-1+I_RANGE+1
 875             Tadv(i,j,k)=t(i,j,k)
 876             if (rmask(i,j).lt.0.5) then
 877               mask_min(i,j)=100.0_r8
 878               Tadv(i,j,k)=1.0
 879               Ta(i,j,k)=1.0
 880             else
 881               mask_min(i,j)=1.0_r8
 882             endif
 883           END DO
 884         END DO
 885       END DO

and skip masked points while calculating Max and Min.

!  Compute UP and DOWN beta-ratios.
 888 !
 889       DO j=J_RANGE
 890         k=1
 891         DO i=I_RANGE
 892           Tmax=MAX(Ta(i-1,j  ,k  )*rmask(i-1,j  ),Tadv(i-1,j  ,k  )     &
 893      &                                              *rmask(i-1,j  ),    &
 894      &             Ta(i  ,j  ,k  )*rmask(i  ,j  ),Tadv(i  ,j  ,k  )     &
 895      &                                              *rmask(i  ,j  ),    &
 896      &             Ta(i+1,j  ,k  )*rmask(i+1,j  ),Tadv(i+1,j  ,k  )     &
 897      &                                              *rmask(i+1,j  ),    &
 898      &             Ta(i  ,j-1,k  )*rmask(i  ,j-1),Tadv(i  ,j-1,k  )     &
 899      &                                              *rmask(i  ,j-1),    &
 900      &             Ta(i  ,j+1,k  )*rmask(i  ,j+1),Tadv(i  ,j+1,k  )     &
 901      &                                              *rmask(i  ,j+1),    &
 902      &             Ta(i  ,j  ,k+1)*rmask(i  ,j  ),Tadv(i  ,j  ,k+1)     &
 903                                                     *rmask(i  ,j  ))
 904           cff1=Ta(i-1,j  ,k  )*MAX(0.0_r8,Ua(i  ,j  ,k  ))-             &
 905      &         Ta(i+1,j  ,k  )*MIN(0.0_r8,Ua(i+1,j  ,k  ))+             &
 906      &         Ta(i  ,j-1,k  )*MAX(0.0_r8,Va(i  ,j  ,k  ))-             &
 907      &         Ta(i  ,j+1,k  )*MIN(0.0_r8,Va(i  ,j+1,k  ))-             &
 908      &         Ta(i  ,j  ,k+1)*MIN(0.0_r8,Wa(i  ,j  ,k  ))
 909           beta_up(i,j,k)=(Tmax-Ta(i,j,k))/(cff1+eps)
 910 !
 911           Tmin=MIN(Ta(i-1,j  ,k  )*mask_min(i-1,j  ),Tadv(i-1,j  ,k  )  &
 912      &                                              *mask_min(i-1,j  ), &
 913      &             Ta(i  ,j  ,k  )*mask_min(i  ,j  ),Tadv(i  ,j  ,k  )  &
 914      &                                              *mask_min(i  ,j  ), &
 915      &             Ta(i+1,j  ,k  )*mask_min(i+1,j  ),Tadv(i+1,j  ,k  )  &
 916      &                                              *mask_min(i+1,j  ), &
 917      &             Ta(i  ,j-1,k  )*mask_min(i  ,j-1),Tadv(i  ,j-1,k  )  &
 918      &                                              *mask_min(i  ,j-1), &
 919      &             Ta(i  ,j+1,k  )*mask_min(i  ,j+1),Tadv(i  ,j+1,k  )  &
 920      &                                              *mask_min(i  ,j+1), &
 921      &             Ta(i  ,j  ,k+1)*mask_min(i  ,j  ),Tadv(i  ,j  ,k+1)  &
 922      &                                              *mask_min(i  ,j  ))
 923           cff2=Ta(i  ,j  ,k  )*MAX(0.0_r8,Ua(i+1,j  ,k  ))-             &
 924      &         Ta(i  ,j  ,k  )*MIN(0.0_r8,Ua(i  ,j  ,k  ))+             &
 925      &         Ta(i  ,j  ,k  )*MAX(0.0_r8,Va(i  ,j+1,k  ))-             &
 926      &         Ta(i  ,j  ,k  )*MIN(0.0_r8,Va(i  ,j  ,k  ))+             &
927      &         Ta(i  ,j  ,k  )*MAX(0.0_r8,Wa(i  ,j  ,k  ))
 928           beta_dn(i,j,k)=(Ta(i,j,k)-Tmin)/(cff2+eps)
 929         END DO
 930         DO k=2,N(ng)-1
 931           DO i=I_RANGE
 932             Tmax=MAX(Ta(i-1,j  ,k  )*rmask(i-1,j  ),Tadv(i-1,j  ,k  )   &
 933                                                     *rmask(i-1,j  ),    &
 934      &               Ta(i  ,j  ,k  )*rmask(i  ,j  ),Tadv(i  ,j  ,k  )   &
 935      &                                              *rmask(i  ,j  ),    &
 936      &               Ta(i+1,j  ,k  )*rmask(i+1,j  ),Tadv(i+1,j  ,k  )   &
 937      &                                              *rmask(i+1,j  ),    &
 938      &               Ta(i  ,j-1,k  )*rmask(i  ,j-1),Tadv(i  ,j-1,k  )   &
 939      &                                              *rmask(i  ,j-1),    &
 940      &               Ta(i  ,j+1,k  )*rmask(i  ,j+1),Tadv(i  ,j+1,k  )   &
 941      &                                              *rmask(i  ,j+1),    &
 942      &               Ta(i  ,j  ,k-1)*rmask(i  ,j  ),Tadv(i  ,j  ,k-1)   &
 943      &                                              *rmask(i  ,j  ),    &
 944      &               Ta(i  ,j  ,k+1)*rmask(i  ,j  ),Tadv(i  ,j  ,k+1)   &
 945      &                                              *rmask(i  ,j  ))
 946             cff1=Ta(i-1,j  ,k  )*MAX(0.0_r8,Ua(i  ,j  ,k  ))-           &
 947      &           Ta(i+1,j  ,k  )*MIN(0.0_r8,Ua(i+1,j  ,k  ))+           &
 948      &           Ta(i  ,j-1,k  )*MAX(0.0_r8,Va(i  ,j  ,k  ))-           &
 949      &           Ta(i  ,j+1,k  )*MIN(0.0_r8,Va(i  ,j+1,k  ))+           &
 950      &           Ta(i  ,j  ,k-1)*MAX(0.0_r8,Wa(i  ,j  ,k-1))-           &
 951      &           Ta(i  ,j  ,k+1)*MIN(0.0_r8,Wa(i  ,j  ,k  ))
 952             beta_up(i,j,k)=(Tmax-Ta(i,j,k))/(cff1+eps)
 953 !
 954             Tmin=MIN(Ta(i-1,j  ,k  )*mask_min(i-1,j  ),Tadv(i-1,j  ,k)  &
 955      &                                              *mask_min(i-1,j  ), &
 956      &               Ta(i  ,j  ,k  )*mask_min(i  ,j  ),Tadv(i  ,j  ,k)  &
 957      &                                              *mask_min(i  ,j  ), &
 958      &               Ta(i+1,j  ,k  )*mask_min(i+1,j  ),Tadv(i+1,j  ,k)  &
 959      &                                              *mask_min(i+1,j  ), &
 960      &               Ta(i  ,j-1,k  )*mask_min(i  ,j-1),Tadv(i  ,j-1,k)  &
 961      &                                              *mask_min(i  ,j-1), &
 962      &               Ta(i  ,j+1,k  )*mask_min(i  ,j+1),Tadv(i  ,j+1,k)  &
 963      &                                              *mask_min(i  ,j+1), &
 964      &               Ta(i  ,j  ,k-1)*mask_min(i  ,j  ),Tadv(i  ,j,k-1)  &
 965      &                                              *mask_min(i  ,j   ),&
 966      &               Ta(i  ,j  ,k+1)*mask_min(i  ,j  ),Tadv(i  ,j,k+1)  &
 967      &                                              *mask_min(i  ,j  ))
 968             cff2=Ta(i  ,j  ,k  )*MAX(0.0_r8,Ua(i+1,j  ,k  ))-           &
 969      &           Ta(i  ,j  ,k  )*MIN(0.0_r8,Ua(i  ,j  ,k  ))+           &
 970      &           Ta(i  ,j  ,k  )*MAX(0.0_r8,Va(i  ,j+1,k  ))-           &
 971      &           Ta(i  ,j  ,k  )*MIN(0.0_r8,Va(i  ,j  ,k  ))+           &
 972      &           Ta(i  ,j  ,k  )*MAX(0.0_r8,Wa(i  ,j  ,k  ))-           &
 973      &           Ta(i  ,j  ,k  )*MIN(0.0_r8,Wa(i  ,j  ,k-1))
 974             beta_dn(i,j,k)=(Ta(i,j,k)-Tmin)/(cff2+eps)
 975           END DO
 976         END DO
 977         k=N(ng)
 978         DO i=I_RANGE
 979           Tmax=MAX(Ta(i-1,j  ,k  )*rmask(i-1,j  ),Tadv(i-1,j  ,k  )     &
 980      &                                              *rmask(i-1,j  ),    &
 981      &             Ta(i  ,j  ,k  )*rmask(i  ,j  ),Tadv(i  ,j  ,k  )     &
 982      &                                              *rmask(i  ,j  ),    &
 983      &             Ta(i+1,j  ,k  )*rmask(i+1,j  ),Tadv(i+1,j  ,k  )     &
 984      &                                              *rmask(i+1,j  ),    &
 985      &             Ta(i  ,j-1,k  )*rmask(i  ,j-1),Tadv(i  ,j-1,k  )     &
 986      &                                              *rmask(i  ,j-1),    &
 987      &             Ta(i  ,j+1,k  )*rmask(i  ,j+1),Tadv(i  ,j+1,k  )     &
 988      &                                              *rmask(i  ,j+1),    &
 989      &             Ta(i  ,j  ,k-1)*rmask(i  ,j  ),Tadv(i  ,j  ,k-1)     &
 990      &                                              *rmask(i  ,j  ))
 991           cff1=Ta(i-1,j  ,k  )*MAX(0.0_r8,Ua(i  ,j  ,k  ))-             &
 992      &         Ta(i+1,j  ,k  )*MIN(0.0_r8,Ua(i+1,j  ,k  ))+             &
 993      &         Ta(i  ,j-1,k  )*MAX(0.0_r8,Va(i  ,j  ,k  ))-             &
 994      &         Ta(i  ,j+1,k  )*MIN(0.0_r8,Va(i  ,j+1,k  ))+             &
 995      &         Ta(i  ,j  ,k-1)*MAX(0.0_r8,Wa(i  ,j  ,k-1))
 996           beta_up(i,j,k)=(Tmax-Ta(i,j,k))/(cff1+eps)
 997 !
 998           Tmin=MIN(Ta(i-1,j  ,k  )*mask_min(i-1,j  ),Tadv(i-1,j  ,k  )  &
 999      &                                              *mask_min(i-1,j  ), &
1000      &             Ta(i  ,j  ,k  )*mask_min(i  ,j  ),Tadv(i  ,j  ,k  )  &
1001      &                                              *mask_min(i  ,j  ), &
1002      &             Ta(i+1,j  ,k  )*mask_min(i+1,j  ),Tadv(i+1,j  ,k  )  &
1003      &                                              *mask_min(i+1,j  ), &
1004      &             Ta(i  ,j-1,k  )*mask_min(i  ,j-1),Tadv(i  ,j-1,k  )  &
1005      &                                              *mask_min(i  ,j-1), &
1006      &             Ta(i  ,j+1,k  )*mask_min(i  ,j+1),Tadv(i  ,j+1,k  )  &
1007      &                                              *mask_min(i  ,j+1), &
1008      &             Ta(i  ,j  ,k-1)*mask_min(i  ,j  ),Tadv(i  ,j  ,k-1)  &
1009      &                                              *mask_min(i  ,j  ))
1010           cff2=Ta(i  ,j  ,k  )*MAX(0.0_r8,Ua(i+1,j  ,k  ))-             &
1011      &         Ta(i  ,j  ,k  )*MIN(0.0_r8,Ua(i  ,j  ,k  ))+             &
1012      &         Ta(i  ,j  ,k  )*MAX(0.0_r8,Va(i  ,j+1,k  ))-             &
1013      &         Ta(i  ,j  ,k  )*MIN(0.0_r8,Va(i  ,j  ,k  ))-             &
1014      &         Ta(i  ,j  ,k  )*MIN(0.0_r8,Wa(i  ,j  ,k-1))
1015           beta_dn(i,j,k)=(Ta(i,j,k)-Tmin)/(cff2+eps)
1016         END DO
1017       END DO
1018 
1019       DO k=1,N(ng)
1020         DO j=J_RANGE
1021           DO i=I_RANGE
1022             if (rmask(i,j).lt.0.5) then
1023                 beta_up(i,j,k)=2.0
1024                 beta_dn(i,j,k)=2.0
1025             endif
1026         END DO
1027        END DO
1028       END DO

Above solution seems to work for my modified GRAV_ADJ case, but I am not a developer and of course I don't know the code very well, I just hide Ta, Tadv, beta_up and beta_dn variables from the Max and Min functions for the masked regions at rho-points.


Adil Sozer (Phd student)
Metu Institute of Marine Sciences, Turkey

[arango]

Interesting, thank you. I have to think about this and make some tests too. The mask_min=100 is problematic for me. We have to find a better way for this. Recall that the tracer in ROMS are not just temperature and salinity. We can have all kind of passive tracers. The code above considers that the MASKING is always on. This code will fail when land/sea masking is not activated. We need to use the MASKING directive here.

[jcwarner]

would your fix work if you only did the mask_min part?
it seems that the masking may not need to be applied to Tadv, since that is set to = t (3), which should already be masked.
I m thinkig to have a mask array that is set to fillvalue for mask=0, and only apply that to the MIN function.need to think on this a little more.
-j

[jcwarner]

Adil-
sorry for the long delay. I am trying to reproduce the problem here, but can not. Perhaps i am not setting the masking in a way that causes the problem. Can you send to me the grav_adj.h file and any ana_mask?? files you modified? I am trying other test cases to see this problem, but i can not re-create it. so i need your help to re-create the problem.
-j
jcwarner@usgs.gov

[adil]

sorry for the late reply,

I havent checked the page for a long time. I will ofcourse help you for reproducing the problem, give me a couple of days in order to find my test cases. I didnt mention before, but the bug is really sensitive to the resolution and disappears with smaller grid size in lock-exchange direction, may be this is why you didnt re-create it.

I will send you an e-mail as soon as possible.


Adil

[adil]

Dr Warner,

I sent my test case files to jcwarner@usgs.gov, may you please inform me when you get my e-mail.

Adil

[ilicakme]

Hi,

Did anybody have a chance to look at this issue? I also have a very simple setup similar
to lock-exchange and I want to use MPDATA as a monotonic scheme, however I'm getting
smaller and larger values than minimum and maximum initial temperatures, respectively.
I'm using WALL options rather than MASKING.

Thanks a lot,
Mehmet

[jcwarner]

yes. i have been poking at this. let me dig deeper.