Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:40:42 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
afe658afe3 Andr*0002        subroutine step_fizhi_corr (myTime, myIter, myThid, dt)
e337e4ca8c Andr*0003 c----------------------------------------------------------------------
                0004 c  Subroutine step_fizhi_corr - 'Wrapper' routine to advance
4420d2948f Jean*0005 c        the physics state and make the new value.
e337e4ca8c Andr*0006 c        At this point, increment with the "correction term"
                0007 c        which includes the dynamics tendency and the integral
                0008 c        constraint to enforce agreement with the dynamics state
                0009 c        Also: Set up "bi, bj loop" and some timers and clocks here.
                0010 c
4420d2948f Jean*0011 c Call:phys2dyn (4) (interpolate physics state to dynamics grid
e337e4ca8c Andr*0012 c                           for use in the correction terms)
                0013 c      AtoC         (convert physics state on dynamics grid to C-Grid)
                0014 c      CtoA         (convert correction term on dynamics grid to A-Grid)
                0015 c      dyn2phys (4) (interpolate A-Grid correction term to physics grid)
                0016 c      step_physics (advance physics state by correction term)
                0017 c-----------------------------------------------------------------------
                0018        implicit none
                0019 #include "SIZE.h"
                0020 #include "GRID.h"
                0021 #include "fizhi_SIZE.h"
f4a0368053 Andr*0022 #include "fizhi_land_SIZE.h"
e337e4ca8c Andr*0023 #include "DYNVARS.h"
                0024 #include "fizhi_coms.h"
                0025 #include "gridalt_mapping.h"
                0026 #include "EEPARAMS.h"
                0027 #include "SURFACE.h"
afe658afe3 Andr*0028 #ifdef ALLOW_DIAGNOSTICS
                0029 #include "fizhi_SHP.h"
                0030 #endif
e337e4ca8c Andr*0031 
3768927683 Andr*0032        integer myIter, myThid
                0033        _RL myTime
afe658afe3 Andr*0034 #ifdef ALLOW_DIAGNOSTICS
                0035       logical  diagnostics_is_on
                0036       external diagnostics_is_on
                0037 #endif
e337e4ca8c Andr*0038 
                0039 c pe on dynamics and physics grid refers to bottom edge
                0040        _RL pephy(1-OLx:sNx+Olx,1-Oly:sNy+Oly,Nrphys+1,nSx,nSy)
                0041        _RL pedyn(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr+1,nSx,nSy)
                0042        _RL windphy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy)
                0043        _RL udyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
                0044        _RL vdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
                0045        _RL thdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
c2377c66e1 Andr*0046        _RL sdyntemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
e337e4ca8c Andr*0047        _RL uphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy)
                0048        _RL vphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy)
                0049        _RL thphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy)
                0050        _RL sphytemp(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy)
d0b11e35fb Andr*0051        _RL tempphy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys,nSx,nSy)
e337e4ca8c Andr*0052 
                0053        integer i, j, L, Lbotij, bi, bj
                0054        integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
f98d2ec0f4 Andr*0055        _RL dt
e337e4ca8c Andr*0056 
5c1a81edc4 Andr*0057        _RL tempij(sNx,sNy)
afe658afe3 Andr*0058        _RL dtfake
732e6393e7 Andr*0059        _RL dtinv
5c1a81edc4 Andr*0060 
e337e4ca8c Andr*0061        im1 = 1-OLx
                0062        im2 = sNx+OLx
                0063        jm1 = 1-OLy
                0064        jm2 = sNy+OLy
                0065        idim1 = 1
                0066        idim2 = sNx
                0067        jdim1 = 1
                0068        jdim2 = sNy
732e6393e7 Andr*0069        dtfake = 1. _d 0
                0070        dtinv = 1. _d 0 / dt
e337e4ca8c Andr*0071 
                0072        do bj = myByLo(myThid), myByHi(myThid)
                0073        do bi = myBxLo(myThid), myBxHi(myThid)
                0074 
                0075 C Construct Pressures on physics and dynamics grids
                0076         do j = 1,sNy
                0077         do i = 1,sNx
                0078          do L = 1,Nr
                0079           pedyn(i,j,L,bi,bj) = 0.
                0080          enddo
                0081         enddo
                0082         enddo
                0083         do j = 1,sNy
                0084         do i = 1,sNx
                0085          Lbotij = ksurfC(i,j,bi,bj)
4420d2948f Jean*0086          if(Lbotij.lt.nr+1)
e337e4ca8c Andr*0087      .    pedyn(i,j,Lbotij,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj)
                0088         enddo
                0089         enddo
                0090         do j = 1,sNy
                0091         do i = 1,sNx
                0092          Lbotij = ksurfC(i,j,bi,bj)
                0093          do L = Lbotij+1,Nr+1
                0094           pedyn(i,j,L,bi,bj) = pedyn(i,j,L-1,bi,bj) -
911d4ccee9 Andr*0095      .               drF(L-1)* rStarExpC(i,j,bi,bj)*hfacC(i,j,L-1,bi,bj)
e337e4ca8c Andr*0096          enddo
                0097 c Do not use a zero field as the top edge pressure for interpolation
                0098          if(pedyn(i,j,Nr+1,bi,bj).lt.1.e-5)
                0099      .                               pedyn(i,j,Nr+1,bi,bj) = 1.e-5
                0100         enddo
                0101         enddo
                0102 
                0103         do j = 1,sNy
                0104         do i = 1,sNx
                0105          pephy(i,j,1,bi,bj)=Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj)
                0106          do L = 2,Nrphys+1
                0107           pephy(i,j,L,bi,bj)=pephy(i,j,L-1,bi,bj)-dpphys(i,j,L-1,bi,bj)
                0108          enddo
                0109 c Do not use a zero field as the top edge pressure for interpolation
                0110          if(pephy(i,j,Nrphys+1,bi,bj).lt.1.e-5)
                0111      .                               pephy(i,j,Nrphys+1,bi,bj) = 1.e-5
                0112         enddo
                0113         enddo
                0114 c
4420d2948f Jean*0115 c Create a wind magnitude field on the physics grid -
732e6393e7 Andr*0116 c   (Load the wind speed bottom up for use by dyn2phys)
e337e4ca8c Andr*0117         do L = 1,Nrphys
                0118         do j = 1,sNy
                0119         do i = 1,sNx
4420d2948f Jean*0120          windphy(i,j,L,bi,bj) =
                0121      .     sqrt(uphy(i,j,Nrphys+1-L,bi,bj)*uphy(i,j,Nrphys+1-L,bi,bj)
732e6393e7 Andr*0122      .        + vphy(i,j,Nrphys+1-L,bi,bj)*vphy(i,j,Nrphys+1-L,bi,bj))
e337e4ca8c Andr*0123         enddo
                0124         enddo
                0125         enddo
180c245802 Andr*0126        enddo
                0127        enddo
                0128 
                0129        CALL TIMER_START('PHYS2DYN          [STEP_FIZHI_CORR]',mythid)
                0130        do bj = myByLo(myThid), myByHi(myThid)
                0131        do bi = myBxLo(myThid), myBxHi(myThid)
e337e4ca8c Andr*0132 
                0133 c Compute correction term (new dyn state-phys state to dyn) on physics grid:
                0134 c    First: interp physics state to dynamics grid
d0b11e35fb Andr*0135 C Note: physics field levels are numbered top down - need bottom up
                0136         do L = 1,Nrphys
                0137         do j = 1,sNy
                0138         do i = 1,sNx
                0139          tempphy(i,j,Nrphys+1-L,bi,bj) = uphy(i,j,L,bi,bj)
                0140         enddo
                0141         enddo
                0142         enddo
                0143         call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,
e337e4ca8c Andr*0144      .        1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,udyntemp)
d0b11e35fb Andr*0145         do L = 1,Nrphys
                0146         do j = 1,sNy
                0147         do i = 1,sNx
                0148          tempphy(i,j,Nrphys+1-L,bi,bj) = vphy(i,j,L,bi,bj)
                0149         enddo
                0150         enddo
                0151         enddo
                0152         call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,
e337e4ca8c Andr*0153      .        1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,vdyntemp)
d0b11e35fb Andr*0154         do L = 1,Nrphys
                0155         do j = 1,sNy
                0156         do i = 1,sNx
                0157          tempphy(i,j,Nrphys+1-L,bi,bj) = thphy(i,j,L,bi,bj)
                0158         enddo
                0159         enddo
                0160         enddo
                0161         call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,
e337e4ca8c Andr*0162      .        1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,thdyntemp)
d0b11e35fb Andr*0163         do L = 1,Nrphys
                0164         do j = 1,sNy
                0165         do i = 1,sNx
                0166          tempphy(i,j,Nrphys+1-L,bi,bj) = sphy(i,j,L,bi,bj)
                0167         enddo
                0168         enddo
                0169         enddo
                0170         call phys2dyn(tempphy,pephy,im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,
e337e4ca8c Andr*0171      .        1,sNx,1,sNy,bi,bj,pedyn,ksurfC,Nr,nlperdyn,sdyntemp)
                0172 
                0173        enddo
                0174        enddo
180c245802 Andr*0175        CALL TIMER_STOP('PHYS2DYN          [STEP_FIZHI_CORR]',mythid)
e337e4ca8c Andr*0176 
                0177 c    Second: Convert physics state on dynamics grid to C-Grid
                0178 
180c245802 Andr*0179        CALL TIMER_START('ATOC              [STEP_FIZHI_CORR]',mythid)
e337e4ca8c Andr*0180        call AtoC(myThid,udyntemp,vdyntemp,maskC,im1,im2,jm1,jm2,Nr,
                0181      .                     Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)
180c245802 Andr*0182        CALL TIMER_STOP('ATOC              [STEP_FIZHI_CORR]',mythid)
e337e4ca8c Andr*0183 
                0184 c    Third: Subtract Phys state on dyn. grid from new dynamics state
                0185        do bj = myByLo(myThid), myByHi(myThid)
                0186        do bi = myBxLo(myThid), myBxHi(myThid)
                0187 
                0188         do L = 1,Nr
                0189         do j = jdim1,jdim2
                0190         do i = idim1,idim2
                0191         udyntemp(i,j,L,bi,bj)=uvel(i,j,L,bi,bj)-udyntemp(i,j,L,bi,bj)
                0192         vdyntemp(i,j,L,bi,bj)=vvel(i,j,L,bi,bj)-vdyntemp(i,j,L,bi,bj)
                0193         thdyntemp(i,j,L,bi,bj)=theta(i,j,L,bi,bj)-thdyntemp(i,j,L,bi,bj)
                0194         sdyntemp(i,j,L,bi,bj)=salt(i,j,L,bi,bj)-sdyntemp(i,j,L,bi,bj)
                0195         enddo
                0196         enddo
                0197         enddo
                0198 
                0199        enddo
                0200        enddo
                0201 
                0202 c    Fourth: Convert correction terms to A-Grid
180c245802 Andr*0203        CALL TIMER_START('CTOA              [STEP_FIZHI_CORR]',mythid)
e337e4ca8c Andr*0204         call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2,
                0205      .     Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)
180c245802 Andr*0206        CALL TIMER_STOP('CTOA              [STEP_FIZHI_CORR]',mythid)
e337e4ca8c Andr*0207 
                0208 c    Fifth: Interpolate correction terms to physics grid
180c245802 Andr*0209        CALL TIMER_START('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
e337e4ca8c Andr*0210        do bj = myByLo(myThid), myByHi(myThid)
                0211        do bi = myBxLo(myThid), myBxHi(myThid)
                0212 
                0213         call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
d0b11e35fb Andr*0214      .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy)
                0215 C Note: correction term is now bottom up - needed in top down arrays
                0216         do L = 1,Nrphys
                0217         do j = 1,sNy
                0218         do i = 1,sNx
                0219          uphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
                0220         enddo
                0221         enddo
                0222         enddo
e337e4ca8c Andr*0223         call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
d0b11e35fb Andr*0224      .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy)
                0225         do L = 1,Nrphys
                0226         do j = 1,sNy
                0227         do i = 1,sNx
                0228          vphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
                0229         enddo
                0230         enddo
                0231         enddo
e337e4ca8c Andr*0232         call dyn2phys(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
d0b11e35fb Andr*0233      .     1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
                0234         do L = 1,Nrphys
                0235         do j = 1,sNy
                0236         do i = 1,sNx
                0237          thphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
                0238         enddo
                0239         enddo
                0240         enddo
e337e4ca8c Andr*0241         call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
d0b11e35fb Andr*0242      .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
                0243         do L = 1,Nrphys
                0244         do j = 1,sNy
                0245         do i = 1,sNx
                0246          sphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
                0247         enddo
                0248         enddo
                0249         enddo
180c245802 Andr*0250        enddo
                0251        enddo
                0252        CALL TIMER_STOP('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
e337e4ca8c Andr*0253 
f98d2ec0f4 Andr*0254 c    Last: Increment physics state by the correction term
180c245802 Andr*0255        do bj = myByLo(myThid), myByHi(myThid)
                0256        do bi = myBxLo(myThid), myBxHi(myThid)
afe658afe3 Andr*0257         call step_physics(uphy,vphy,thphy,sphy,dtfake,im1,im2,jm1,jm2,
f98d2ec0f4 Andr*0258      .   Nrphys,Nsx,Nsy,1,sNx,1,sNy,bi,bj,
                0259      .                            uphytemp,vphytemp,thphytemp,sphytemp)
e337e4ca8c Andr*0260 
8598abfc35 Andr*0261         call qcheck (im1,im2,jm1,jm2,Nrphys,Nsx,Nsy,idim1,idim2,
                0262      .         jdim1,jdim2,bi,bj,dpphys,sphy)
                0263 
                0264        enddo
                0265        enddo
                0266 
732e6393e7 Andr*0267 #ifdef ALLOW_DIAGNOSTICS
61fe9f1bad Andr*0268       call diagnostics_fill(uphy,'UAVE    ',0,Nrphys,0,1,1,myThid)
e5d2509d4f Andr*0269       call diagnostics_fill(vphy,'VAVE    ',0,Nrphys,0,1,1,myThid)
61fe9f1bad Andr*0270       call diagnostics_fill(thphy,'TAVE    ',0,Nrphys,0,1,1,myThid)
                0271       call diagnostics_fill(sphy,'QAVE    ',0,Nrphys,0,1,1,myThid)
                0272 #endif
                0273 
                0274 #ifdef ALLOW_DIAGNOSTICS
afe658afe3 Andr*0275       do bj = myByLo(myThid), myByHi(myThid)
                0276       do bi = myBxLo(myThid), myBxHi(myThid)
                0277       do L=1,Nrphys
8598abfc35 Andr*0278 
afe658afe3 Andr*0279 c Total Tendency on Fizhi Grid for U (m/sec/day)
                0280 c -----------------------------------------------
                0281        if(diagnostics_is_on('TENDUFIZ',myThid) ) then
                0282         do j=jm1,jm2
                0283         do i=im1,im2
                0284          tempij(i,j) = (uphy (i,j,L,bi,bj)-ubef(i,j,L,bi,bj) )
732e6393e7 Andr*0285      .                     * 86400. _d 0 * dtinv
afe658afe3 Andr*0286         enddo
                0287         enddo
                0288         call diagnostics_fill(tempij,'TENDUFIZ',L,1,3,bi,bj,myThid)
                0289        endif
                0290 
                0291 c Total Tendency on Fizhi Grid for V (m/sec/day)
                0292 c -----------------------------------------------
                0293        if(diagnostics_is_on('TENDVFIZ',myThid) ) then
                0294         do j=jm1,jm2
                0295         do i=im1,im2
                0296          tempij(i,j) = (vphy (i,j,L,bi,bj)-vbef(i,j,L,bi,bj) )
732e6393e7 Andr*0297      .                     * 86400. _d 0 * dtinv
afe658afe3 Andr*0298         enddo
                0299         enddo
                0300         call diagnostics_fill(tempij,'TENDVFIZ',L,1,3,bi,bj,myThid)
                0301        endif
                0302 
                0303 c Total Tendency on Fizhi Grid for U (m/sec/day)
                0304 c -----------------------------------------------
                0305        if(diagnostics_is_on('TENDTFIZ',myThid) ) then
                0306         do j=jm1,jm2
                0307         do i=im1,im2
                0308          tempij(i,j) = (thphy (i,j,L,bi,bj)-thbef(i,j,L,bi,bj) )
732e6393e7 Andr*0309      .                     * 86400. _d 0 * dtinv
afe658afe3 Andr*0310         enddo
                0311         enddo
                0312         call diagnostics_fill(tempij,'TENDTFIZ',L,1,3,bi,bj,myThid)
                0313        endif
                0314 
                0315 c Total Tendency on Fizhi Grid for U (m/sec/day)
                0316 c -----------------------------------------------
                0317        if(diagnostics_is_on('TENDQFIZ',myThid) ) then
                0318         do j=jm1,jm2
                0319         do i=im1,im2
                0320          tempij(i,j) = (sphy (i,j,L,bi,bj)-sbef(i,j,L,bi,bj) )
732e6393e7 Andr*0321      .                     * 86400. _d 0 * dtinv
afe658afe3 Andr*0322         enddo
                0323         enddo
                0324         call diagnostics_fill(tempij,'TENDQFIZ',L,1,3,bi,bj,myThid)
                0325        endif
                0326 
                0327       enddo
                0328       enddo
                0329       enddo
e337e4ca8c Andr*0330 
732e6393e7 Andr*0331 c Gridalt Correction Term Tendency for U and V (m/sec/day)
                0332 c --------------------------------------------------------
                0333       if(diagnostics_is_on('CORRDU  ',myThid) .or.
                0334      .    diagnostics_is_on('CORRDV  ',myThid)   ) then
                0335 
                0336 C gridalt correction term - first step is to compute adv+filters tendency
                0337 C                           on dynamics grid (total - physics tend)
                0338       do bj = myByLo(myThid), myByHi(myThid)
                0339       do bi = myBxLo(myThid), myBxHi(myThid)
                0340        do L=1,Nr
                0341        do j=jm1,jm2
                0342        do i=im1,im2
4420d2948f Jean*0343         udyntemp(i,j,L,bi,bj) =
                0344      .    (uvel(i,j,L,bi,bj)-udynbef(i,j,L,bi,bj))*dtinv -
732e6393e7 Andr*0345      .                                        guphy(i,j,L,bi,bj)
4420d2948f Jean*0346         vdyntemp(i,j,L,bi,bj) =
                0347      .    (vvel(i,j,L,bi,bj)-vdynbef(i,j,L,bi,bj))*dtinv -
732e6393e7 Andr*0348      .                                        gvphy(i,j,L,bi,bj)
                0349        enddo
                0350        enddo
                0351        enddo
                0352 C Next step - interpolate to fizhi grid
                0353 C  first put the u and v tendencies on an a-grid
                0354        CALL TIMER_START('CTOA              [STEP_FIZHI_CORR]',mythid)
                0355        call CtoA(myThid,udyntemp,vdyntemp,maskW,maskS,im1,im2,jm1,jm2,
                0356      .     Nr,Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)
                0357        CALL TIMER_STOP('CTOA              [STEP_FIZHI_CORR]',mythid)
                0358 C  then do vertical interpolation
                0359        do L = 1,Nrphys
                0360        do j = 1,sNy
                0361        do i = 1,sNx
4420d2948f Jean*0362         windphy(i,j,L,bi,bj) =
                0363      .     sqrt(uphy(i,j,Nrphys+1-L,bi,bj)*uphy(i,j,Nrphys+1-L,bi,bj)
732e6393e7 Andr*0364      .        + vphy(i,j,Nrphys+1-L,bi,bj)*vphy(i,j,Nrphys+1-L,bi,bj))
                0365        enddo
                0366        enddo
                0367        enddo
                0368        CALL TIMER_START('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
                0369         call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
                0370      .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy)
                0371        CALL TIMER_STOP('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
                0372 C Note: adv+filters term is now bottom up - needed in top down arrays
                0373        do L = 1,Nrphys
                0374        do j = 1,sNy
                0375        do i = 1,sNx
                0376          uphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
                0377        enddo
                0378        enddo
                0379        enddo
                0380        call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
                0381      .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy)
                0382        do L = 1,Nrphys
                0383        do j = 1,sNy
                0384        do i = 1,sNx
                0385          vphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
                0386        enddo
                0387        enddo
                0388        enddo
                0389 C Last Step - subtract adv+filters and physics tend from total tend on physics grid
                0390        do L = 1,Nrphys
                0391        do j = 1,sNy
                0392        do i = 1,sNx
                0393         uphytemp(i,j,L,bi,bj)=
                0394      .    ( (uphy(i,j,L,bi,bj)-ubef(i,j,L,bi,bj))*dtinv
                0395      .     - duphy(i,j,L,bi,bj) - uphytemp(i,j,L,bi,bj) ) * 86400. _d 0
                0396         vphytemp(i,j,L,bi,bj)=
                0397      .    ( (vphy(i,j,L,bi,bj)-vbef(i,j,L,bi,bj))*dtinv
                0398      .     - dvphy(i,j,L,bi,bj) - vphytemp(i,j,L,bi,bj) ) * 86400. _d 0
                0399        enddo
                0400        enddo
                0401        enddo
                0402       enddo
                0403       enddo
                0404 
                0405       if(diagnostics_is_on('CORRDU  ',myThid)) then
                0406        call diagnostics_fill(uphytemp,'CORRDU  ',0,Nrphys,0,1,1,myThid)
                0407       endif
                0408       if(diagnostics_is_on('CORRDV  ',myThid)) then
                0409        call diagnostics_fill(vphytemp,'CORRDV  ',0,Nrphys,0,1,1,myThid)
                0410       endif
                0411 
                0412       endif
                0413 
                0414 c Gridalt Correction Term Tendency for TH (deg K/day)
                0415 c --------------------------------------------------------
                0416       if(diagnostics_is_on('CORRDT  ',myThid))  then
                0417 
                0418 C gridalt correction term - first step is to compute adv+filters tendency
                0419 C                           on dynamics grid (total - physics tend)
                0420       do bj = myByLo(myThid), myByHi(myThid)
                0421       do bi = myBxLo(myThid), myBxHi(myThid)
                0422        do L=1,Nr
                0423        do j=jm1,jm2
                0424        do i=im1,im2
4420d2948f Jean*0425         thdyntemp(i,j,L,bi,bj) =
                0426      .    (theta(i,j,L,bi,bj)-thdynbef(i,j,L,bi,bj))*dtinv -
732e6393e7 Andr*0427      .                                        gthphy(i,j,L,bi,bj)
                0428        enddo
                0429        enddo
                0430        enddo
                0431 C Next step - interpolate to fizhi grid
                0432        CALL TIMER_START('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
                0433        call dyn2phys(thdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
                0434      .     1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
                0435        CALL TIMER_STOP('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
                0436 C Note: adv+filters term is now bottom up - needed in top down arrays
                0437        do L = 1,Nrphys
                0438        do j = 1,sNy
                0439        do i = 1,sNx
                0440          thphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
                0441        enddo
                0442        enddo
                0443        enddo
                0444 C Last Step - subtract adv+filters and physics tend from total tend on physics grid
                0445        do L = 1,Nrphys
                0446        do j = 1,sNy
                0447        do i = 1,sNx
                0448         thphytemp(i,j,L,bi,bj)=
                0449      .    ( (thphy(i,j,L,bi,bj)-thbef(i,j,L,bi,bj))*dtinv
4420d2948f Jean*0450      .     - dthphy(i,j,L,bi,bj) - thphytemp(i,j,L,bi,bj)
                0451      .    ) * 86400. _d 0
732e6393e7 Andr*0452        enddo
                0453        enddo
                0454        enddo
                0455       enddo
                0456       enddo
                0457 
                0458       call diagnostics_fill(thphytemp,'CORRDT  ',0,Nrphys,0,1,1,myThid)
                0459       endif
                0460 
                0461 c Gridalt Correction Term Tendency for Q (kg/kg/day)
                0462 c --------------------------------------------------------
                0463       if(diagnostics_is_on('CORRDQ  ',myThid))  then
                0464 
                0465 C gridalt correction term - first step is to compute adv+filters tendency
                0466 C                           on dynamics grid (total - physics tend)
                0467       do bj = myByLo(myThid), myByHi(myThid)
                0468       do bi = myBxLo(myThid), myBxHi(myThid)
                0469        do L=1,Nr
                0470        do j=jm1,jm2
                0471        do i=im1,im2
4420d2948f Jean*0472         sdyntemp(i,j,L,bi,bj) =
                0473      .    (salt(i,j,L,bi,bj)-sdynbef(i,j,L,bi,bj))*dtinv -
732e6393e7 Andr*0474      .                                        gsphy(i,j,L,bi,bj)
                0475        enddo
                0476        enddo
                0477        enddo
                0478 C Next step - interpolate to fizhi grid
                0479        CALL TIMER_START('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
                0480        call dyn2phys(sdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,1,sNx,
                0481      .      1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
                0482        CALL TIMER_STOP('DYN2PHYS          [STEP_FIZHI_CORR]',mythid)
                0483 C Note: adv+filters term is now bottom up - needed in top down arrays
                0484        do L = 1,Nrphys
                0485        do j = 1,sNy
                0486        do i = 1,sNx
                0487          sphytemp(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
                0488        enddo
                0489        enddo
                0490        enddo
                0491 C Last Step - subtract adv+filters and physics tend from total tend on physics grid
                0492        do L = 1,Nrphys
                0493        do j = 1,sNy
                0494        do i = 1,sNx
                0495         sphytemp(i,j,L,bi,bj)=
                0496      .    ( (sphy(i,j,L,bi,bj)-sbef(i,j,L,bi,bj))*dtinv
                0497      .     - dsphy(i,j,L,bi,bj) - sphytemp(i,j,L,bi,bj) ) * 86400. _d 0
                0498        enddo
                0499        enddo
                0500        enddo
                0501       enddo
                0502       enddo
                0503 
                0504       call diagnostics_fill(sphytemp,'CORRDQ  ',0,Nrphys,0,1,1,myThid)
                0505       endif
                0506 #endif
                0507 
afe658afe3 Andr*0508       return
                0509       end
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!