140 |
INTEGER nipass,ipass |
INTEGER nipass,ipass |
141 |
INTEGER myTile, nCFace |
INTEGER myTile, nCFace |
142 |
LOGICAL N_edge, S_edge, E_edge, W_edge |
LOGICAL N_edge, S_edge, E_edge, W_edge |
143 |
|
#ifdef ALLOW_DIAGNOSTICS |
144 |
|
CHARACTER*8 diagName |
145 |
|
CHARACTER*4 GAD_DIAG_SUFX, diagSufx |
146 |
|
EXTERNAL GAD_DIAG_SUFX |
147 |
|
#endif |
148 |
CEOP |
CEOP |
149 |
|
|
150 |
#ifdef ALLOW_AUTODIFF_TAMC |
#ifdef ALLOW_AUTODIFF_TAMC |
168 |
endif |
endif |
169 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
170 |
|
|
171 |
|
#ifdef ALLOW_DIAGNOSTICS |
172 |
|
C-- Set diagnostic suffix for the current tracer |
173 |
|
IF ( useDiagnostics ) THEN |
174 |
|
diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid ) |
175 |
|
ENDIF |
176 |
|
#endif |
177 |
|
|
178 |
C-- Set up work arrays with valid (i.e. not NaN) values |
C-- Set up work arrays with valid (i.e. not NaN) values |
179 |
C These inital values do not alter the numerical results. They |
C These inital values do not alter the numerical results. They |
180 |
C just ensure that all memory references are to valid floating |
C just ensure that all memory references are to valid floating |
256 |
ENDDO |
ENDDO |
257 |
ENDDO |
ENDDO |
258 |
|
|
259 |
|
#ifndef ALLOW_AUTODIFF_TAMC |
260 |
IF (useCubedSphereExchange) THEN |
IF (useCubedSphereExchange) THEN |
261 |
withSigns = .FALSE. |
withSigns = .FALSE. |
262 |
CALL FILL_CS_CORNER_UV_RS( |
CALL FILL_CS_CORNER_UV_RS( |
263 |
& withSigns, maskLocW,maskLocS, bi,bj, myThid ) |
& withSigns, maskLocW,maskLocS, bi,bj, myThid ) |
264 |
ENDIF |
ENDIF |
265 |
|
#endif |
266 |
|
|
267 |
C-- Multiple passes for different directions on different tiles |
C-- Multiple passes for different directions on different tiles |
268 |
C-- For cube need one pass for each of red, green and blue axes. |
C-- For cube need one pass for each of red, green and blue axes. |
317 |
C and b) the overlap of myTile are not cube-face Edges |
C and b) the overlap of myTile are not cube-face Edges |
318 |
IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN |
IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN |
319 |
|
|
320 |
|
#ifndef ALLOW_AUTODIFF_TAMC |
321 |
C- Internal exchange for calculations in X |
C- Internal exchange for calculations in X |
322 |
#ifdef MULTIDIM_OLD_VERSION |
#ifdef MULTIDIM_OLD_VERSION |
323 |
IF ( useCubedSphereExchange ) THEN |
IF ( useCubedSphereExchange ) THEN |
327 |
#endif |
#endif |
328 |
CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid ) |
CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid ) |
329 |
ENDIF |
ENDIF |
330 |
|
#endif |
331 |
|
|
332 |
C- Advective flux in X |
C- Advective flux in X |
333 |
DO j=1-Oly,sNy+Oly |
DO j=1-Oly,sNy+Oly |
344 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
345 |
|
|
346 |
IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN |
IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN |
347 |
CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, deltaTtracer, |
CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k), |
348 |
I uTrans, uVel, maskLocW, localTij, |
I uTrans, uVel, maskLocW, localTij, |
349 |
O af, myThid ) |
O af, myThid ) |
350 |
ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN |
ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN |
351 |
CALL GAD_DST3_ADV_X( bi,bj,k, deltaTtracer, |
CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k), |
352 |
I uTrans, uVel, maskLocW, localTij, |
I uTrans, uVel, maskLocW, localTij, |
353 |
O af, myThid ) |
O af, myThid ) |
354 |
ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN |
ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN |
355 |
CALL GAD_DST3FL_ADV_X( bi,bj,k, deltaTtracer, |
CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k), |
356 |
I uTrans, uVel, maskLocW, localTij, |
I uTrans, uVel, maskLocW, localTij, |
357 |
O af, myThid ) |
O af, myThid ) |
358 |
ELSE |
ELSE |
362 |
C- Advective flux in X : done |
C- Advective flux in X : done |
363 |
ENDIF |
ENDIF |
364 |
|
|
365 |
|
#ifndef ALLOW_AUTODIFF_TAMC |
366 |
C- Internal exchange for next calculations in Y |
C- Internal exchange for next calculations in Y |
367 |
IF ( overlapOnly .AND. ipass.EQ.1 ) THEN |
IF ( overlapOnly .AND. ipass.EQ.1 ) THEN |
368 |
CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid ) |
CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid ) |
369 |
ENDIF |
ENDIF |
370 |
|
#endif |
371 |
|
|
372 |
C- Update the local tracer field where needed: |
C- Update the local tracer field where needed: |
373 |
|
|
383 |
IF ( S_edge ) THEN |
IF ( S_edge ) THEN |
384 |
DO j=1-Oly,0 |
DO j=1-Oly,0 |
385 |
DO i=iMinUpd,iMaxUpd |
DO i=iMinUpd,iMaxUpd |
386 |
localTij(i,j)=localTij(i,j)-deltaTtracer* |
localTij(i,j)=localTij(i,j)-dTtracerLev(k)* |
387 |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
388 |
& *recip_rA(i,j,bi,bj) |
& *recip_rA(i,j,bi,bj) |
389 |
& *( af(i+1,j)-af(i,j) |
& *( af(i+1,j)-af(i,j) |
395 |
IF ( N_edge ) THEN |
IF ( N_edge ) THEN |
396 |
DO j=sNy+1,sNy+Oly |
DO j=sNy+1,sNy+Oly |
397 |
DO i=iMinUpd,iMaxUpd |
DO i=iMinUpd,iMaxUpd |
398 |
localTij(i,j)=localTij(i,j)-deltaTtracer* |
localTij(i,j)=localTij(i,j)-dTtracerLev(k)* |
399 |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
400 |
& *recip_rA(i,j,bi,bj) |
& *recip_rA(i,j,bi,bj) |
401 |
& *( af(i+1,j)-af(i,j) |
& *( af(i+1,j)-af(i,j) |
413 |
IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy |
IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy |
414 |
DO j=jMinUpd,jMaxUpd |
DO j=jMinUpd,jMaxUpd |
415 |
DO i=1-Olx+1,sNx+Olx-1 |
DO i=1-Olx+1,sNx+Olx-1 |
416 |
localTij(i,j)=localTij(i,j)-deltaTtracer* |
localTij(i,j)=localTij(i,j)-dTtracerLev(k)* |
417 |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
418 |
& *recip_rA(i,j,bi,bj) |
& *recip_rA(i,j,bi,bj) |
419 |
& *( af(i+1,j)-af(i,j) |
& *( af(i+1,j)-af(i,j) |
435 |
CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid ) |
CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid ) |
436 |
ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN |
ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN |
437 |
CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid ) |
CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid ) |
438 |
|
#ifdef ALLOW_PTRACERS |
439 |
|
ELSEIF (tracerIdentity.GE.GAD_TR1) THEN |
440 |
|
CALL OBCS_APPLY_PTRACER( bi, bj, k, |
441 |
|
& tracerIdentity-GAD_TR1+1, localTij, myThid ) |
442 |
|
#endif /* ALLOW_PTRACERS */ |
443 |
ENDIF |
ENDIF |
444 |
ENDIF |
ENDIF |
445 |
#endif /* ALLOW_OBCS */ |
#endif /* ALLOW_OBCS */ |
459 |
C and b) the overlap of myTile are not cube-face edges |
C and b) the overlap of myTile are not cube-face edges |
460 |
IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN |
IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN |
461 |
|
|
462 |
|
#ifndef ALLOW_AUTODIFF_TAMC |
463 |
C- Internal exchange for calculations in Y |
C- Internal exchange for calculations in Y |
464 |
#ifdef MULTIDIM_OLD_VERSION |
#ifdef MULTIDIM_OLD_VERSION |
465 |
IF ( useCubedSphereExchange ) THEN |
IF ( useCubedSphereExchange ) THEN |
469 |
#endif |
#endif |
470 |
CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid ) |
CALL FILL_CS_CORNER_TR_RL(.FALSE., localTij, bi,bj, myThid ) |
471 |
ENDIF |
ENDIF |
472 |
|
#endif |
473 |
|
|
474 |
C- Advective flux in Y |
C- Advective flux in Y |
475 |
DO j=1-Oly,sNy+Oly |
DO j=1-Oly,sNy+Oly |
486 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
487 |
|
|
488 |
IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN |
IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN |
489 |
CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, deltaTtracer, |
CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k), |
490 |
I vTrans, vVel, maskLocS, localTij, |
I vTrans, vVel, maskLocS, localTij, |
491 |
O af, myThid ) |
O af, myThid ) |
492 |
ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN |
ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN |
493 |
CALL GAD_DST3_ADV_Y( bi,bj,k, deltaTtracer, |
CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k), |
494 |
I vTrans, vVel, maskLocS, localTij, |
I vTrans, vVel, maskLocS, localTij, |
495 |
O af, myThid ) |
O af, myThid ) |
496 |
ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN |
ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN |
497 |
CALL GAD_DST3FL_ADV_Y( bi,bj,k, deltaTtracer, |
CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k), |
498 |
I vTrans, vVel, maskLocS, localTij, |
I vTrans, vVel, maskLocS, localTij, |
499 |
O af, myThid ) |
O af, myThid ) |
500 |
ELSE |
ELSE |
504 |
C- Advective flux in Y : done |
C- Advective flux in Y : done |
505 |
ENDIF |
ENDIF |
506 |
|
|
507 |
|
#ifndef ALLOW_AUTODIFF_TAMC |
508 |
C- Internal exchange for next calculations in X |
C- Internal exchange for next calculations in X |
509 |
IF ( overlapOnly .AND. ipass.EQ.1 ) THEN |
IF ( overlapOnly .AND. ipass.EQ.1 ) THEN |
510 |
CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid ) |
CALL FILL_CS_CORNER_TR_RL( .TRUE., localTij, bi,bj, myThid ) |
511 |
ENDIF |
ENDIF |
512 |
|
#endif |
513 |
|
|
514 |
C- Update the local tracer field where needed: |
C- Update the local tracer field where needed: |
515 |
|
|
525 |
IF ( W_edge ) THEN |
IF ( W_edge ) THEN |
526 |
DO j=jMinUpd,jMaxUpd |
DO j=jMinUpd,jMaxUpd |
527 |
DO i=1-Olx,0 |
DO i=1-Olx,0 |
528 |
localTij(i,j)=localTij(i,j)-deltaTtracer* |
localTij(i,j)=localTij(i,j)-dTtracerLev(k)* |
529 |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
530 |
& *recip_rA(i,j,bi,bj) |
& *recip_rA(i,j,bi,bj) |
531 |
& *( af(i,j+1)-af(i,j) |
& *( af(i,j+1)-af(i,j) |
537 |
IF ( E_edge ) THEN |
IF ( E_edge ) THEN |
538 |
DO j=jMinUpd,jMaxUpd |
DO j=jMinUpd,jMaxUpd |
539 |
DO i=sNx+1,sNx+Olx |
DO i=sNx+1,sNx+Olx |
540 |
localTij(i,j)=localTij(i,j)-deltaTtracer* |
localTij(i,j)=localTij(i,j)-dTtracerLev(k)* |
541 |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
542 |
& *recip_rA(i,j,bi,bj) |
& *recip_rA(i,j,bi,bj) |
543 |
& *( af(i,j+1)-af(i,j) |
& *( af(i,j+1)-af(i,j) |
555 |
IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx |
IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx |
556 |
DO j=1-Oly+1,sNy+Oly-1 |
DO j=1-Oly+1,sNy+Oly-1 |
557 |
DO i=iMinUpd,iMaxUpd |
DO i=iMinUpd,iMaxUpd |
558 |
localTij(i,j)=localTij(i,j)-deltaTtracer* |
localTij(i,j)=localTij(i,j)-dTtracerLev(k)* |
559 |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
560 |
& *recip_rA(i,j,bi,bj) |
& *recip_rA(i,j,bi,bj) |
561 |
& *( af(i,j+1)-af(i,j) |
& *( af(i,j+1)-af(i,j) |
577 |
CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid ) |
CALL OBCS_APPLY_TLOC( bi, bj, k, localTij, myThid ) |
578 |
ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN |
ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN |
579 |
CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid ) |
CALL OBCS_APPLY_SLOC( bi, bj, k, localTij, myThid ) |
580 |
|
#ifdef ALLOW_PTRACERS |
581 |
|
ELSEIF (tracerIdentity.GE.GAD_TR1) THEN |
582 |
|
CALL OBCS_APPLY_PTRACER( bi, bj, k, |
583 |
|
& tracerIdentity-GAD_TR1+1, localTij, myThid ) |
584 |
|
#endif /* ALLOW_PTRACERS */ |
585 |
ENDIF |
ENDIF |
586 |
ENDIF |
ENDIF |
587 |
#endif /* ALLOW_OBCS */ |
#endif /* ALLOW_OBCS */ |
600 |
DO j=1-Oly,sNy+Oly |
DO j=1-Oly,sNy+Oly |
601 |
DO i=1-Olx,sNx+Olx |
DO i=1-Olx,sNx+Olx |
602 |
gTracer(i,j,k,bi,bj)= |
gTracer(i,j,k,bi,bj)= |
603 |
& (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer |
& (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k) |
604 |
ENDDO |
ENDDO |
605 |
ENDDO |
ENDDO |
606 |
ELSE |
ELSE |
612 |
ENDDO |
ENDDO |
613 |
ENDIF |
ENDIF |
614 |
|
|
615 |
|
#ifdef ALLOW_DIAGNOSTICS |
616 |
|
IF ( useDiagnostics ) THEN |
617 |
|
diagName = 'ADVx'//diagSufx |
618 |
|
CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid) |
619 |
|
diagName = 'ADVy'//diagSufx |
620 |
|
CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid) |
621 |
|
ENDIF |
622 |
|
#endif |
623 |
|
|
624 |
#ifdef ALLOW_DEBUG |
#ifdef ALLOW_DEBUG |
625 |
IF ( debugLevel .GE. debLevB |
IF ( debugLevel .GE. debLevB |
626 |
& .AND. tracerIdentity.EQ.GAD_TEMPERATURE |
& .AND. tracerIdentity.EQ.GAD_TEMPERATURE |
699 |
C- Compute vertical advective flux in the interior: |
C- Compute vertical advective flux in the interior: |
700 |
IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN |
IF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN |
701 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
702 |
CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTtracer, |
CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, dTtracerLev(k), |
703 |
I rTrans, wVel, localTijk, |
I rTrans, wVel, localTijk, |
704 |
O fVerT(1-Olx,1-Oly,kUp), myThid ) |
O fVerT(1-Olx,1-Oly,kUp), myThid ) |
705 |
ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN |
ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN |
706 |
CALL GAD_DST3_ADV_R( bi,bj,k, deltaTtracer, |
CALL GAD_DST3_ADV_R( bi,bj,k, dTtracerLev(k), |
707 |
I rTrans, wVel, localTijk, |
I rTrans, wVel, localTijk, |
708 |
O fVerT(1-Olx,1-Oly,kUp), myThid ) |
O fVerT(1-Olx,1-Oly,kUp), myThid ) |
709 |
ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN |
ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN |
710 |
CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTtracer, |
CALL GAD_DST3FL_ADV_R( bi,bj,k, dTtracerLev(k), |
711 |
I rTrans, wVel, localTijk, |
I rTrans, wVel, localTijk, |
712 |
O fVerT(1-Olx,1-Oly,kUp), myThid ) |
O fVerT(1-Olx,1-Oly,kUp), myThid ) |
713 |
ELSE |
ELSE |
727 |
C-- Divergence of vertical fluxes |
C-- Divergence of vertical fluxes |
728 |
DO j=1-Oly,sNy+Oly |
DO j=1-Oly,sNy+Oly |
729 |
DO i=1-Olx,sNx+Olx |
DO i=1-Olx,sNx+Olx |
730 |
localTij(i,j)=localTijk(i,j,k)-deltaTtracer* |
localTij(i,j)=localTijk(i,j,k)-dTtracerLev(k)* |
731 |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
& _recip_hFacC(i,j,k,bi,bj)*recip_drF(k) |
732 |
& *recip_rA(i,j,bi,bj) |
& *recip_rA(i,j,bi,bj) |
733 |
& *( fVerT(i,j,kUp)-fVerT(i,j,kDown) |
& *( fVerT(i,j,kDown)-fVerT(i,j,kUp) |
734 |
& -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j)) |
& -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j)) |
735 |
& )*rkFac |
& )*rkSign |
736 |
gTracer(i,j,k,bi,bj)= |
gTracer(i,j,k,bi,bj)= |
737 |
& (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer |
& (localTij(i,j)-tracer(i,j,k,bi,bj))/dTtracerLev(k) |
738 |
ENDDO |
ENDDO |
739 |
ENDDO |
ENDDO |
740 |
|
|
741 |
|
#ifdef ALLOW_DIAGNOSTICS |
742 |
|
IF ( useDiagnostics ) THEN |
743 |
|
diagName = 'ADVr'//diagSufx |
744 |
|
CALL DIAGNOSTICS_FILL( fVerT(1-Olx,1-Oly,kUp), |
745 |
|
& diagName, k,1, 2,bi,bj, myThid) |
746 |
|
ENDIF |
747 |
|
#endif |
748 |
|
|
749 |
C-- End of K loop for vertical flux |
C-- End of K loop for vertical flux |
750 |
ENDDO |
ENDDO |
751 |
C-- end of if not.implicitAdvection block |
C-- end of if not.implicitAdvection block |