C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/shelfice_remeshing/code/Attic/shelfice_remeshing.F,v 1.1 2015/07/21 14:52:40 dgoldberg Exp $ C $Name: $ #include "SHELFICE_OPTIONS.h" #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif #ifdef ALLOW_CTRL # include "CTRL_OPTIONS.h" #endif CBOP C !ROUTINE: SHELFICE_THERMODYNAMICS C !INTERFACE: SUBROUTINE SHELFICE_REMESHING( I myTime, myIter, myThid ) C !DESCRIPTION: \bv C *=============================================================* C | S/R SHELFICE_THERMODYNAMICS C | o shelf-ice main routine. C | compute temperature and (virtual) salt flux at the C | shelf-ice ocean interface C | C | stresses at the ice/water interface are computed in separate C | routines that are called from mom_fluxform/mom_vecinv C *=============================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "DYNVARS.h" #include "FFIELDS.h" #include "SHELFICE.h" #include "SHELFICE_COST.h" #ifdef ALLOW_AUTODIFF # include "CTRL_SIZE.h" # include "ctrl.h" # include "ctrl_dummy.h" #endif /* ALLOW_AUTODIFF */ #ifdef ALLOW_AUTODIFF_TAMC # ifdef SHI_ALLOW_GAMMAFRICT # include "tamc.h" # include "tamc_keys.h" # endif /* SHI_ALLOW_GAMMAFRICT */ #endif /* ALLOW_AUTODIFF_TAMC */ C === Functions ==== LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myIter :: iteration counter for this thread C myTime :: time counter for this thread C myThid :: thread number for this instance of the routine. _RL myTime INTEGER myIter INTEGER myThid _RL hFacCtmp _RL hFacMnSz _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C === Local variables === C I,J,K,Kp1,bi,bj :: loop counters C tLoc, sLoc, pLoc :: local in-situ temperature, salinity, pressure C theta/saltFreeze :: temperature and salinity of water at the C ice-ocean interface (at the freezing point) C freshWaterFlux :: local variable for fresh water melt flux due C to melting in kg/m^2/s C (negative density x melt rate) C convertFW2SaltLoc:: local copy of convertFW2Salt C cFac :: 1 for conservative form, 0, otherwise C rFac :: realFreshWaterFlux factor C dFac :: 0 for diffusive heat flux (Holland and Jenkins, 1999, C eq21) C 1 for advective and diffusive heat flux (eq22, 26, 31) C fwflxFac :: only effective for dFac=1, 1 if we expect a melting C fresh water flux, 0 otherwise C auxiliary variables and abbreviations: C a0, a1, a2, b, c0 C eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8 C aqe, bqe, cqe, discrim, recip_aqe C drKp1, recip_drLoc INTEGER I,J,K,Kp1 INTEGER bi,bj _RL tLoc(1:sNx,1:sNy) _RL sLoc(1:sNx,1:sNy) _RL pLoc(1:sNx,1:sNy) _RL uLoc(1:sNx,1:sNy) _RL vLoc(1:sNx,1:sNy) _RL u_topdr(1:sNx+1,1:sNy+1,nSx,nSy) _RL v_topdr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL thetaFreeze, saltFreeze, recip_Cp _RL freshWaterFlux, convertFW2SaltLoc _RL a0, a1, a2, b, c0 _RL eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8 _RL cFac, rFac, dFac, fwflxFac _RL aqe, bqe, cqe, discrim, recip_aqe _RL drKp1, recip_drLoc _RL recip_latentHeat _RL tmpFac _RL shiPr, shiSc, shiLo, recip_shiKarman, shiTwoThirds _RL gammaTmoleT, gammaTmoleS, gammaTurb, gammaTurbConst _RL ustar, ustarSq, etastar CEOP C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| IF ( SHELFICERemeshFrequency .EQ. myTime) THEN print *, 'JJ WOZ HERE',R_shelfIce DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO J = 1, sNy+1 DO I = 1, sNx+1 IF (etaN(I,J,bi,bj) .GT. 4.4 ) THEN K = MAX(1,kTopC(I,J,bi,bj)) c salt(I,J,bi,bj,K-1)=salt(I,J,bi,bj,K) c theta(I,J,bi,bj,K-1)=theta(I,J,bi,bj,K) c uVel(I,J,bi,bj,K-1)=uVel(I,J,bi,bj,K) c uVel(I+1,J,bi,bj,K-1)=uVel(I+1,J,bi,bj,K) c vVel(I,J,bi,bj,K-1)=uVel(I,J,bi,bj,K) c vVel(I,J+1,bi,bj,K-1)=uVel(I,J+1,bi,bj,K) etaN(I,J,bi,bj) = etaN(I,J,bi,bj) - 10.0 R_shelfIce(I,J,bi,bj) = R_shelfIce(I,J,bi,bj) + 10.0 ENDIF ENDDO ENDDO ENDDO ENDDO DO K=1, Nr hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) ) DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx c o Non-dimensional distance between grid boundary and model surface hFacCtmp = (rF(k)-R_shelfIce(I,J,bi,bj))*recip_drF(K) C o Reduce the previous fraction : substract the outside part. hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0) c o set to zero if empty Column : hFacCtmp = max( hFacCtmp, 0. _d 0) C o Impose minimum fraction and/or size (dimensional) IF (hFacCtmp.LT.hFacMnSz) THEN IF (hFacCtmp.LT.hFacMnSz*0.5) THEN hFacC(I,J,K,bi,bj)=0. ELSE hFacC(I,J,K,bi,bj)=hFacMnSz ENDIF ELSE hFacC(I,J,K,bi,bj)=hFacCtmp ENDIF ENDDO ENDDO ENDDO c cC- Re-calculate Reference surface position, taking into account hFacC cC initialize Total column fluid thickness and surface k index cC Note: if no fluid (continent) ==> kSurf = Nr+1 DO bj=myByLo(myThid), myByHi(myThid) DO bi=myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx tmpfld(i,j,bi,bj) = 0. kSurfC(i,j,bi,bj) = Nr+1 c maskH(i,j,bi,bj) = 0. Ro_surf(i,j,bi,bj) = R_low(i,j,bi,bj) DO k=Nr,1,-1 Ro_surf(i,j,bi,bj)=Ro_surf(i,j,bi,bj) & +drF(k)*hFacC(i,j,k,bi,bj) IF (hFacC(i,j,k,bi,bj).NE.0.) THEN kSurfC(i,j,bi,bj) = k c maskH(i,j,bi,bj) = 1. tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1. ENDIF ENDDO kLowC(i,j,bi,bj) = 0 DO k= 1, Nr IF (hFacC(i,j,k,bi,bj).NE.0) THEN kLowC(i,j,bi,bj) = k ENDIF ENDDO maskInC(i,j,bi,bj)= 0. C Weird IF loop here JJ IF ( kSurfC(i,j,bi,bj).LE.Nr ) THEN maskInC(i,j,bi,bj)= 1. ENDIF ENDDO ENDDO C- end bi,bj loops. ENDDO ENDDO c c IF ( printDomain ) THEN cc CALL PLOT_FIELD_XYRS( tmpfld, c & 'Model Depths K Index' , -1, myThid ) c CALL PLOT_FIELD_XYRS(R_low, c & 'Model R_low (ini_masks_etc)', -1, myThid ) c CALL PLOT_FIELD_XYRS(Ro_surf, c & 'Model Ro_surf (ini_masks_etc)', -1, myThid ) c ENDIF c cC-- Calculate quantities derived from XY depth map DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx C Total fluid column thickness (r_unit) : c Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) C Inverse of fluid column thickness (1/r_unit) IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN recip_Rcol(i,j,bi,bj) = 0. ELSE recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpfld(i,j,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDDO c cC-- hFacW and hFacS (at U and V points) DO bj=myByLo(myThid), myByHi(myThid) DO bi=myBxLo(myThid), myBxHi(myThid) DO k=1, Nr DO j=1-OLy,sNy+OLy hFacW(1-OLx,j,k,bi,bj)= 0. DO i=2-OLx,sNx+OLx hFacW(i,j,k,bi,bj)= & MIN(hFacC(i,j,k,bi,bj),hFacC(i-1,j,k,bi,bj)) ENDDO ENDDO DO i=1-OLx,sNx+OLx hFacS(i,1-OLy,k,bi,bj)= 0. ENDDO DO j=2-OLy,sNy+oly DO i=1-OLx,sNx+OLx hFacS(i,j,k,bi,bj)= & MIN(hFacC(i,j,k,bi,bj),hFacC(i,j-1,k,bi,bj)) ENDDO ENDDO ENDDO c cC rLow & reference rSurf at Western & Southern edges (U and V points) i = 1-OLx DO j=1-OLy,sNy+OLy rLowW (i,j,bi,bj) = 0. rSurfW(i,j,bi,bj) = 0. ENDDO j = 1-OLy DO i=1-OLx,sNx+OLx rLowS (i,j,bi,bj) = 0. rSurfS(i,j,bi,bj) = 0. ENDDO DO j=1-OLy,sNy+OLy DO i=2-OLx,sNx+OLx rLowW(i,j,bi,bj) = & MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) ) rSurfW(i,j,bi,bj) = & MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) ) rSurfW(i,j,bi,bj) = & MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) ) ENDDO ENDDO DO j=2-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rLowS(i,j,bi,bj) = & MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) ) rSurfS(i,j,bi,bj) = & MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) ) rSurfS(i,j,bi,bj) = & MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) ) ENDDO ENDDO C- end bi,bj loops. ENDDO ENDDO CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid) CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid ) CALL EXCH_UV_XY_RS( rLowW, rLowS, .FALSE., myThid ) c cC-- Addtional closing of Western and Southern grid-cell edges: for example, cC a) might add some "thin walls" in specific location cC-- b) close non-periodic N & S boundaries of lat-lon grid at the N/S poles. CALL ADD_WALLS2MASKS( myThid ) cC-- Calculate surface k index for interface W & S (U & V points) DO bj=myByLo(myThid), myByHi(myThid) DO bi=myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx kSurfW(i,j,bi,bj) = Nr+1 kSurfS(i,j,bi,bj) = Nr+1 DO k=Nr,1,-1 IF (hFacW(i,j,k,bi,bj).NE.0.) kSurfW(i,j,bi,bj) = k IF (hFacS(i,j,k,bi,bj).NE.0.) kSurfS(i,j,bi,bj) = k ENDDO maskInW(i,j,bi,bj)= 0. cC Wierd if statements JJ c IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN maskInW(i,j,bi,bj)= 1. maskInS(i,j,bi,bj)= 0. ENDIF IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN maskInS(i,j,bi,bj)= 1. ENDIF ENDDO ENDDO ENDDO ENDDO c ELSE #ifndef DISABLE_SIGMA_CODE C--- Sigma and Hybrid-Sigma set-up: CALL INI_SIGMA_HFAC( myThid ) c#endif /* DISABLE_SIGMA_CODE */ cc ENDIF cC---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Write to disk: Total Column Thickness & hFac(C,W,S): C This I/O is now done in write_grid.F c CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid) c CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid) c CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid) c CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid) c IF ( printDomain ) THEN CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 0, myThid ) CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 0, myThid ) CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 0, myThid ) ENDIF C-- Masks and reciprocals of hFac[CWS] DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF (hFacC(i,j,k,bi,bj) .NE. 0. ) THEN c recip_hFacC(i,j,k,bi,bj) = 1. _d 0 / hFacC(i,j,k,bi,bj) c maskC(i,j,k,bi,bj) = 1. ELSE c recip_hFacC(i,j,k,bi,bj) = 0. c maskC(i,j,k,bi,bj) = 0. ENDIF IF (hFacW(i,j,k,bi,bj) .NE. 0. ) THEN c recip_hFacW(i,j,k,bi,bj) = 1. _d 0 / hFacW(i,j,k,bi,bj) c maskW(i,j,k,bi,bj) = 1. ELSE c recip_hFacW(i,j,k,bi,bj) = 0. c maskW(i,j,k,bi,bj) = 0. ENDIF IF (hFacS(i,j,k,bi,bj) .NE. 0. ) THEN c recip_hFacS(i,j,k,bi,bj) = 1. _d 0 / hFacS(i,j,k,bi,bj) c maskS(i,j,k,bi,bj) = 1. ELSE c recip_hFacS(i,j,k,bi,bj) = 0. c maskS(i,j,k,bi,bj) = 0. ENDIF ENDDO ENDDO ENDDO c#ifdef NONLIN_FRSURF C-- Save initial geometrical hFac factor into h0Fac (fixed in time): C Note: In case 1 pkg modifies hFac (from packages_init_fixed, called C later in sequence of calls) this pkg would need also to update h0Fac. c DO k=1,Nr c DO j=1-OLy,sNy+OLy c DO i=1-OLx,sNx+OLx c h0FacC(i,j,k,bi,bj) = _hFacC(i,j,k,bi,bj) c h0FacW(i,j,k,bi,bj) = _hFacW(i,j,k,bi,bj) c h0FacS(i,j,k,bi,bj) = _hFacS(i,j,k,bi,bj) c ENDDO c ENDDO c ENDDO #endif /* NONLIN_FRSURF */ C- end bi,bj loops. ENDDO ENDDO ENDIF RETURN END