Bugs in pkg/layer: Does pkg/layers work in pressure coordinates? and more.

[mjlosch]

• Comparing pkg/layers output for meridional overturning circulation computations between pressure and z-coordinates, I have the impression that pkg/layers code does not work properly in pressure coordinates. In the code, there are references to w-points above and below, directions from "warm" to "cold", which has the opposite meaning in pressure coordinates than in z-coordinates and to the surface level as k=1, e.g.

MITgcm/pkg/layers/layers_thermodynamics.F

Lines 67 to 87 in f1a2968

	k = 1
	C -- Loop for surface fluxes
	DO j=1-OLy,sNy+OLy
	DO i=1-OLx,sNx+OLx
	#ifdef SHORTWAVE_HEATING
	C -- Have to remove the shortwave from the surface flux because it is added later
	IF (iTracer.EQ.1) THEN
	layers_surfflux(i,j,k,iTracer,bi,bj) =
	& layers_surfflux(i,j,k,iTracer,bi,bj)
	C -- Sign convention for Qsw means we have to add it to subtract it
	& +Qsw(i,j,bi,bj)
	ENDIF
	#endif /* SHORTWAVE_HEATING */
	layers_surfflux(i,j,k,iTracer,bi,bj) =
	& layers_surfflux(i,j,k,iTracer,bi,bj)
	& *recip_drF(1)*_recip_hFacC(i,j,1,bi,bj)
	& *fluxfac(iTracer)
	ENDDO
	ENDDO

or

MITgcm/pkg/layers/layers_fluxcalc.F

Lines 476 to 485 in f1a2968

	C Initialize the search indices
	DO j = 1,sNy+1
	DO i = 1,sNx+1
	C The temperature index (layer_G) goes from cold to warm.
	C The water column goes from warm (k=1) to cold (k=Nr).
	C So initialize the search with the warmest value.
	kgc(i,j) = Nlayers
	kgcw(i,j) = Nlayers-1
	ENDDO
	ENDDO

or

MITgcm/pkg/layers/layers_fluxcalc.F

Lines 582 to 586 in f1a2968

	C find the layer indices kgcw for the w point
	CALL LAYERS_LOCATE(
	I layers_bounds_w(1,iLa),Nlayers-1,mSteps,sNx,sNy,TatC,
	O kgcw,
	I myThid )

• Further, I have the impression that the following lines of code are not doing what they are supposed to do (replacing the .EQV. operator, which is most likely my own fault):

MITgcm/pkg/layers/layers_fluxcalc.F

Lines 765 to 767 in f1a2968

	CML IF ((xx(n).GE.xx(1)).EQV.(x(i,j).GE.xx(km))) THEN
	IF ( ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))).OR.
	& ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))) ) THEN

and should probably be

           IF ( ( xx(n).GE.xx(1).AND.x(i,j).GE.xx(km) ).OR. 
      &         ( xx(n).LT.xx(1).AND.x(i,j).LT.xx(km) ) ) THEN 

[jm-c]

@mjlosch I have 2 comments:

1. There might be some problems with layers_thermodynamics.F (e.g., issue Layers package budget #60) but this piece of code is not used to compute the meridional overturning in density coordinate (routine is empty when LAYERS_THERMODYNAMICS is "#undef"). So, in some sense, it's a different issue.
2. The layers_fluxcalc.F has been used to get the transport in potential density coordinate or in temperature coordinate, with different ordering direction (temp increasing from bottom to top but density decreasing from bottom to top), but some comments might not be accurate. Also, I vaguely remember (but could be wrong) that it has also been used for the atmosphere in P-coordinate. This does not mean the problem you saw is not related to p-coord but its not that the implementation is completely missing.