C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/dgoldberg/streamice/streamice_advect_thickness.F,v 1.10 2013/06/12 20:48:08 dgoldberg Exp $ C $Name: $ #include "STREAMICE_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP SUBROUTINE STREAMICE_ADVECT_THICKNESS ( myThid,myIter,time_step ) C /============================================================\ C | SUBROUTINE | C | o | C |============================================================| C | | C \============================================================/ IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "GRID.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "STREAMICE.h" #include "STREAMICE_ADV.h" #ifdef ALLOW_AUTODIFF_TAMC # include "tamc.h" #endif INTEGER myThid, myIter _RL time_step #ifdef ALLOW_STREAMICE INTEGER i, j, bi, bj, Gi, Gj _RL thick_bd, uflux, vflux, max_icfl, loc_icfl _RL time_step_full, time_step_rem _RL sec_per_year, time_step_loc, MR, SMB, TMB _RL BCVALX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BCVALY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS BCMASKX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS BCMASKY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL utrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vtrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL h_after_uflux_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL h_after_vflux_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL hflux_x_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL hflux_y_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CHARACTER*(MAX_LEN_MBUF) msgBuf sec_per_year = 365.*86400. time_step_loc = time_step / sec_per_year time_step_full = time_step_loc time_step_rem = time_step_loc PRINT *, "time_step_loc ", time_step_loc #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics #endif DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy+1 DO i=1,sNx+1 H_streamice_prev(i,j,bi,bj) = & H_streamice(i,j,bi,bj) hflux_x_SI (i,j,bi,bj) = 0. _d 0 hflux_y_SI (i,j,bi,bj) = 0. _d 0 IF (STREAMICE_ufacemask(i,j,bi,bj).eq.3.0) THEN BCMASKX(i,j,bi,bj) = 3.0 BCVALX(i,j,bi,bj) = h_ubdry_values_SI(i,j,bi,bj) utrans(i,j,bi,bj) = .5 * ( & u_streamice(i,j,bi,bj)+u_streamice(i,j+1,bi,bj)) ELSEIF (STREAMICE_ufacemask(i,j,bi,bj).eq.4.0) THEN uflux = u_flux_bdry_SI(i,j,bi,bj) BCMASKX(i,j,bi,bj) = 0.0 BCVALX(i,j,bi,bj) = uflux utrans(i,j,bi,bj) = 1.0 ELSEIF (.not.( & STREAMICE_hmask(i,j,bi,bj).eq.1.0.OR. & STREAMICE_hmask(i-1,j,bi,bj).eq.1.0)) THEN BCMASKX(i,j,bi,bj) = 0.0 BCVALX(i,j,bi,bj) = 0. _d 0 utrans(i,j,bi,bj) = 0. _d 0 ELSE BCMASKX(i,j,bi,bj) = 0.0 BCVALX(i,j,bi,bj) = 0. _d 0 utrans(i,j,bi,bj) = .5 * ( & u_streamice(i,j,bi,bj)+u_streamice(i,j+1,bi,bj)) ENDIF IF (STREAMICE_vfacemask(i,j,bi,bj).eq.3.0) THEN BCMASKy(i,j,bi,bj) = 3.0 BCVALy(i,j,bi,bj) = h_vbdry_values_SI(i,j,bi,bj) vtrans(i,j,bi,bj) = .5 * ( & v_streamice(i,j,bi,bj)+v_streamice(i+1,j,bi,bj)) ELSEIF (STREAMICE_vfacemask(i,j,bi,bj).eq.4.0) THEN vflux = v_flux_bdry_SI(i,j,bi,bj) BCMASKy(i,j,bi,bj) = 0.0 BCVALy(i,j,bi,bj) = vflux vtrans(i,j,bi,bj) = 1.0 ELSEIF (.not.( & STREAMICE_hmask(i,j,bi,bj).eq.1.0.OR. & STREAMICE_hmask(i,j-1,bi,bj).eq.1.0)) THEN BCMASKY(i,j,bi,bj) = 0.0 BCVALY(i,j,bi,bj) = 0. _d 0 vtrans(i,j,bi,bj) = 0. _d 0 ELSE BCMASKy(i,j,bi,bj) = 0.0 BCVALy(i,j,bi,bj) = 0. _d 0 vtrans(i,j,bi,bj) = .5 * ( & v_streamice(i,j,bi,bj)+v_streamice(i+1,j,bi,bj)) ENDIF ENDDO ENDDO ENDDO ENDDO _EXCH_XY_RL(utrans,myThid) _EXCH_XY_RL(vtrans,myThid) _EXCH_XY_RS(BCMASKx,myThid) _EXCH_XY_RS(BCMASKy,myThid) _EXCH_XY_RL(BCVALX,myThid) _EXCH_XY_RL(BCVALY,myThid) #ifndef ALLOW_AUTODIFF_TAMC max_icfl = 1.e-20 DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN loc_icfl=max(abs(utrans(i,j,bi,bj)), & abs(utrans(i+1,j,bi,bj))) / dxF(i,j,bi,bj) loc_icfl=max(loc_icfl,max(abs(vtrans(i,j,bi,bj)), & abs(vtrans(i,j+1,bi,bj))) / dyF(i,j,bi,bj)) if (loc_icfl.gt.max_icfl) then max_icfl = loc_icfl ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO CALL GLOBAL_MAX_R8 (max_icfl, myThid) #endif #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics CADJ STORE H_streamice = comlev1, key=ikey_dynamics #endif #ifndef ALLOW_AUTODIFF_TAMC do while (time_step_rem .gt. 1.e-15) time_step_loc = min ( & streamice_cfl_factor / max_icfl, & time_step_rem ) if (time_step_loc .lt. time_step_full) then PRINT *, "TAKING PARTIAL TIME STEP", time_step_loc endif #endif CALL STREAMICE_ADV_FLUX_FL_X ( myThid , I utrans , I H_streamice , I BCMASKX, I BCVALX, O hflux_x_SI, I time_step_loc ) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-3,sNy+3 DO i=1,sNx Gi = (myXGlobalLo-1)+(bi-1)*sNx+i Gj = (myYGlobalLo-1)+(bj-1)*sNy+j IF (((Gj .ge. 1) .and. (Gj .le. Ny)) & .or.STREAMICE_NS_PERIODIC) THEN IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN h_after_uflux_SI (i,j,bi,bj) = H_streamice(i,j,bi,bj) - & (hflux_x_SI(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) - & hflux_x_SI(i,j,bi,bj)*dyG(i,j,bi,bj)) & * recip_rA (i,j,bi,bj) * time_step_loc IF ( h_after_uflux_SI (i,j,bi,bj).le.0.0) THEN PRINT *, "h neg after x", i,j,hflux_x_SI(i+1,j,bi,bj), & hflux_x_SI(i,j,bi,bj) ENDIF ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics #endif ! CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid, ! O hflux_y_SI, ! O h_after_vflux_SI, ! I time_step_loc ) CALL STREAMICE_ADV_FLUX_FL_Y ( myThid , I vtrans , I h_after_uflux_si , I BCMASKY, I BCVALY, O hflux_y_SI, I time_step_loc ) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx Gi = (myXGlobalLo-1)+(bi-1)*sNx+i Gj = (myYGlobalLo-1)+(bj-1)*sNy+j IF (((Gj .ge. 1) .and. (Gj .le. Ny)) & .or.STREAMICE_EW_PERIODIC) THEN IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN h_after_vflux_SI (i,j,bi,bj) = h_after_uflux_SI(i,j,bi,bj) - & (hflux_y_SI(i,j+1,bi,bj)*dxG(i,j+1,bi,bj) - & hflux_y_SI(i,j,bi,bj)*dxG(i,j,bi,bj)) * & recip_rA (i,j,bi,bj) * time_step_loc IF ( h_after_vflux_SI (i,j,bi,bj).le.0.0) THEN PRINT *, "h neg after y", i,j,hflux_y_SI(i,j+1,bi,bj), & hflux_y_SI(i,j,bi,bj) ENDIF ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN H_streamice (i,j,bi,bj) = & h_after_vflux_SI (i,j,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDDO ! NOTE: AT THIS POINT H IS NOT VALID ON OVERLAP!!! ! CALL STREAMICE_ADV_FRONT ( ! & myThid, time_step_loc, ! & hflux_x_SI, hflux_y_SI ) #ifdef ALLOW_STREAMICE_2DTRACER CALL STREAMICE_ADVECT_2DTRACER( & myThid, & myIter, & time_step, & uTrans, & vTrans, & bcMaskx, & bcMasky ) #endif #ifndef ALLOW_AUTODIFF_TAMC time_step_rem = time_step_rem - time_step_loc enddo #endif ! NOW WE APPLY MELT RATES !! ! THIS MAY BE MOVED TO A SEPARATE SUBROUTINE DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0 .or. & STREAMICE_hmask(i,j,bi,bj).eq.2.0) THEN MR = (1.-float_frac_streamice(i,j,bi,bj)) * & BDOT_STREAMICE(i,j,bi,bj) SMB = ADOT_STREAMICE(i,j,bi,bj) TMB = SMB - MR IF ((TMB.lt.0.0) .and. & (MR * time_step_loc .gt. & H_streamice (i,j,bi,bj))) THEN H_streamice (i,j,bi,bj) = 0. _d 0 STREAMICE_hmask(i,j,bi,bj) = 0. ELSE H_streamice (i,j,bi,bj) = & H_streamice (i,j,bi,bj) + TMB * time_step_loc ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO _EXCH_XY_RL (H_streamice,myThid) WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT , 1) #endif END