do l=1,lm

      integer  i,j,l,im1,ip1,imax,m

      do l = 1,Lm

      do l = 2,Lm

       do l = 2,nbase(i)

      do l = 2,nbase(i)

        do l = 1,nthin(i)

       do l = nthin(i)+1,nbase(i)

        if( l.eq.2 ) crif2 = min(0.7 _d 0,crif2)

       do l = nbase(i)+1,Lm+1

        if( l.eq.Lm+1 ) then

      do l = nthin(i)+1,Lm

        do l = 2,nthin(i)

      do l = 1,Lm

      integer i,j,l,n,nn

       rains(i,l) = 0.

      integer i,j,iloop,indx,indgath,l,nn,num,numdeps,nt

        IF((qqz(i,j,l).gt.qqcolmax(i,j))

         qqcolmax(i,j) = qqz(i,j,l)

        if((qqz(i,j,l).lt.qqcolmin(i,j))

      do l=1,lm

      dumoist(i,j,l) = 0.

      dvmoist(i,j,l) = 0.

      dtmoist(i,j,l) = 0.

        dqmoist(i,j,l,nt) = 0.

        CVTH(I,L) =  (TH (I,L) - saveth(i,l))

        CVQ (I,L) =  (SHL(I,L) - saveq (i,l))

        TH (I,L) = saveth(i,l) + cvth(i,l)

        SHL(I,L) = saveq (i,l) + cvq (i,l)

         CVU(I,L,nt) = ( UL(I,L,nt)-saveu(i,l,nt) )*sp(i)*tminv

       TMP1(I,L) = sp(i) * ((TL(I,L)/PLK(I,L))-saveth(i,l)) * tminv

        TMP1(I,L) = ((TL(I,L)/PLK(I,L))-saveth(i,l))*plk(i,l)*sday*tminv

          do l=lup,lm

            if (saveth(i,l)*plk(i,l).le. tice ) then

     . ( lwdt(i,j,l,bi,bj) + 

     .                      * pk(i,j,l,bi,bj)*pinv(i,j)*86400

     . ( lwdtclr(i,j,l,bi,bj) + 

     .                      * pk(i,j,l,bi,bj)*pinv(i,j)*86400

     .  + swdt(i,j,l,bi,bj)*radswt(i,j,bi,bj)*

     .                   pk(i,j,l,bi,bj)*pinv(i,j)*86400

     .  + swdtclr(i,j,l,bi,bj)*radswt(i,j,bi,bj)*

     .                   pk(i,j,l,bi,bj)*pinv(i,j)*86400

     .             moistt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)

     .             turbt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)

     .             *pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)

     .             swdt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)

              alf = min( max((tzl(i,l)-253.15)/20.,0. _d 0) ,1. _d 0)

      do l=1,lm

      do l=1,nlay

      sttflux(i,l) = sttflux(i,l) * cp

      do l=1,nlay

      stqflux(i,l) = stqflux(i,l) * alhl

      do l = 1,nlay

      do l = 1,nlay

          tmpdiag(i,1) = dtturb(i,1,l) * pkz(i,1,l)*atimstp*secday

      do l =1,nlay

        duturb(i,1,l) = duturb(i,1,l)*atimstp

        dvturb(i,1,l) = dvturb(i,1,l)*atimstp

        dtturb(i,1,l) = dtturb(i,1,l)*thtgz(i)

         dqturb(i,1,l,nt) = dqturb(i,1,l,nt)*thtgz(i)

         IF ( (qq(i,l).gt.qqcolmax(I)).and.(levpbl(i).eq.nlev))then

          qqcolmax(i) = qq(i,l)

         if((qq(i,l).lt.qqcolmin(i)).and.(levpbl(i).eq.nlev))

     1                                                levpbl(i)=l

      do  l=1,nlay

      call STRIP2TILE(uz(1,1,l),igrd,tmpstrip,nchp,ijall,

      tends(i) = ( u(i,l)-tmpstrip(i) )

      call PASTIT (tends,dumsc(1,l),istrip,nchp,nchptot,1,nn)

      call STRIP2TILE(vz(1,1,l),igrd,tmpstrip,nchp,ijall,

      tends(i) = ( v(i,l)-tmpstrip(i) )

      call PASTIT (tends,dvmsc(1,l),istrip,nchp,nchptot,1,nn)

      call STRIP2TILE(tz(1,1,l),igrd,tmpstrip,nchp,ijall,

      tends(i) = ( th(i,l)-tmpstrip(i) )

      call PASTIT (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)

      call STRIP2TILE(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,

      tends(i) = ( sh(i,l)-tmpstrip(i) )

      call PASTIT (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)

      subroutine STRIPIT(a,b,irun,ia,ib,l,k)

      integer ia,ib,irun,l,k

      _RL a(ia,l), b(ib,l)

        do 100 j=1,l

        do 200 j=1,l

      subroutine STRIPITINT(a,b,irun,ia,ib,l,k)

      integer ia,ib,irun,l,k,a(ia,l),b(ib,l)

        do 100 j=1,l

        do 200 j=1,l

      _RL a(ia,l), b(ib,l)

      subroutine PSTBITINT(b,a,ib,ia,irun,l,k)

      _RL a(ia,l)

      integer b(ib,l)

      subroutine PSTBMPIT(b,a,ib,ia,irun,l,k)

      _RL a(ia,l), b(ib,l)

       pkz(i,j,L) = ( ple (i,j,l+1)*pkle(i,j,l+1)

     .             - ple (i,j,l)*pkle(i,j,l) )

     .           /  ( akap1* (ple (i,j,l+1)-ple (i,j,l)) )

       print *,' thphy at level ',l,' ',tempij

       print *,' sphy at level ',l,' ',tempij

       print *,' sphy after qcheck at level ',l,' ',tempij

      INTEGER k,l

      integer i,j,k,l, n,n1, ngrid, ds,de

          l = mnc_fd_ind(indf,k)

          n1 = ILNBLNK(mnc_d_names(l))

     &         .AND. (mnc_d_names(l)(1:n1) .EQ. dname(1:n))) THEN

      integer k,l,m,q,bi,bj,iG,jG, myThid,num_file,length_for_tile

         DO l=0,sNx+1

        xy_buffer_r8(l,m)=xC(l,m,bi,bj)

         DO l=0,sNx+1

                xy_buffer_r8(l,m)=yC(l,m,bi,bj)

         do l=1,1000

           profiles_data_buff(m,l,k,bi,bj)=0

           profiles_weight_buff(m,l,k,bi,bj)=0

      integer l, ilinsysinfo, ii, jj

            do l = 1, n

               phtmpadmtlm(l) = workd(ipntr(1)+l-1)

            do l = 1, n

               workd(ipntr(2)+l-1) = phtmpadmtlm(l)

               do l = 1, n

                  tmpvec1(l) = v(l,j)

               call DCOPY( m, tmpvec2, 1, u(l,j), 1 )

               temp = one / dnrm2( m, u(l,j), 1 )

               print *, 'ph-print ', j, dnrm2( m, u(l,j), 1 )

               call DSCAL(m , temp, u(l,j), 1 )

      INTEGER iG,jG,lG, l

      DO l=1,length

         lG= (jG-1)*nPx*length + (iG-1)*length + l

         global(lG) = local(l)

         DO l=1,length

            lG= (jG-1)*nPx*length + (iG-1)*length + l

            global(lG) = local(l)

            DO l=1,length

               lG= (jG-1)*nPx*length + (iG-1)*length + l

               global(lG) = temp(l)

      INTEGER iG,jG,lG, l

      DO l=1,length

         lG= (jG-1)*nPx*length + (iG-1)*length + l

         local(l) = global(lG)

         DO l=1,length

            lG= (jG-1)*nPx*length + (iG-1)*length + l

            local(l) = global(lG)

            DO l=1,length

               lG= (jG-1)*nPx*length + (iG-1)*length + l

               temp(l) = global(lG)

           DO l=1,diagSt_nbLists

            DO k=1,diagSt_nbActv(l)

             IF ( iSdiag(k,l).EQ.-iSp ) THEN

               region2fill(j) = MAX(region2fill(j),diagSt_region(j,l))

      INTEGER iSp, ndId, j,l

           DO l=1,diagSt_nbLists

            DO k=1,diagSt_nbActv(l)

             IF ( iSdiag(k,l).EQ.-iSp ) THEN

               region2fill(j) = MAX(region2fill(j),diagSt_region(j,l))

      INTEGER m, n, j, k, l, bi, bj

      INTEGER i, j, l

       l = MIN(k,sizF)

     &                         *frcFldRL(i,j,l,bi,bj)

     &                         *frcFldRS(i,j,l,bi,bj)

     &                         *frcFldRL(i,j,l,bi,bj)

     &                         *frcFldRS(i,j,l,bi,bj)

       l = MIN(k,sizF)

     &                         *frcFldRL(i,j,l,bi,bj)

     &                         *frcFldRS(i,j,l,bi,bj)

     &                         *frcFldRL(i,j,l,bi,bj)

     &                         *frcFldRS(i,j,l,bi,bj)

           DO l=1,diagSt_nbLists

            DO k=1,diagSt_nbActv(l)

             IF ( iSdiag(k,l).EQ.-iSp ) THEN

               region2fill(j) = MAX(region2fill(j),diagSt_region(j,l))

      INTEGER m, n, j, k, l, bi, bj

           DO l=1,diagSt_nbLists

            DO k=1,diagSt_nbActv(l)

             IF ( iSdiag(k,l).EQ.-iSp ) THEN

               region2fill(j) = MAX(region2fill(j),diagSt_region(j,l))

      INTEGER m, n, j, k, l, bi, bj

      INTEGER  i,l

      DO l = 1,ldimLoc

        frequency(l)  = 0.

        timePhase(l)  = UNSET_RL

        averagingFreq(l) = 0.

        averagingPhase(l)= 0.

        repeatCycle(l)   = 0

        fileName(l)   = blkFilName

        missing_value(l)     = UNSET_RL

        missing_value_int(l) = UNSET_I

        fileFlags(l)  = blk8c

          levels(k,l) = UNSET_RL

          fields(m,l) = blk8c

      DO l = 1,ldimLoc

        stat_freq(l)  = 0.

        stat_phase(l) = UNSET_RL

        stat_fname(l) = blkFilName

          stat_region(k,l) = UNSET_I

          stat_fields(m,l) = blk8c

      DO l = 1,ldimLoc

       iLen = ILNBLNK(fileName(l))

         freq(n)    = frequency(l)

         IF ( timePhase(l).NE. UNSET_RL ) THEN

           phase(n) = timePhase(l)

         ELSEIF ( frequency(l) .LT. 0. ) THEN

           phase(n) = -0.5 _d 0 * frequency(l)

         IF ( averagingFreq(l).GT.0. .AND. repeatCycle(l).GT.1 ) THEN

           averageFreq(n)  = averagingFreq(l)

           averagePhase(n) = averagingPhase(l)

           averageCycle(n) = repeatCycle(l)

         ELSEIF (averagingFreq(l).NE.0. .OR. repeatCycle(l).NE.0) THEN

     &       averagingFreq(l), repeatCycle(l)

     &         ' for list l=', l, ', fileName: ', fileName(l)

         ELSEIF ( frequency(l) .EQ. 0. ) THEN

         ELSEIF ( frequency(l) .GT. 0. ) THEN

           averageFreq(n)  = frequency(l)

         IF ( missing_value(l) .NE. UNSET_RL )

     &        misvalFlt(n) = missing_value(l)

         IF ( missing_value_int(l) .NE. UNSET_I )

     &        misvalInt(n) = missing_value_int(l)

         fnames(n)  = fileName (l)

         fflags(n)  = fileFlags(l)

         IF ( levels(1,l).NE.UNSET_RL ) THEN

             IF ( levels(k,l).NE.UNSET_RL .AND.

               levs(nlevels(n),n) = levels(k,l)

             ELSEIF ( levels(k,l).NE.UNSET_RL ) THEN

     &         'when trying to add level(k=', k, ' )=', levels(k,l)

     &         ' for list l=', l, ', fileName: ', fileName(l)

           IF ( fields(m,l).NE.blk8c .AND.

             flds(nfields(n),n) = fields(m,l)

           ELSEIF ( fields(m,l).NE.blk8c ) THEN

     &        'when trying to add field (m=', m, ' ): ',fields(m,l)

     &        ' in list l=', l, ', fileName: ', fileName(l)

     &    'when trying to add list l=', l

     &    ' Frq=', frequency(l), ', fileName: ', fileName(l)

      DO l = 1,ldimLoc

       iLen = ILNBLNK(stat_fname(l))

         diagSt_freq(n) = stat_freq(l)

         IF ( stat_phase(l).NE. UNSET_RL ) THEN

           diagSt_phase(n) = stat_phase(l)

         ELSEIF ( stat_freq(l) .LT. 0. ) THEN

           diagSt_phase(n) = -0.5 _d 0 * stat_freq(l)

         diagSt_Fname(n)  = stat_fname(l)

           j = stat_region(k,l)

     &        ' in list l=', l, ', stat_fname: ', stat_fname(l)

     &         ' in list l=', l, ', stat_fname: ', stat_fname(l)

           IF ( stat_fields(m,l).NE.blk8c .AND.

             diagSt_Flds(diagSt_nbFlds(n),n) = stat_fields(m,l)

           ELSEIF ( stat_fields(m,l).NE.blk8c ) THEN

     &        ' ): ',stat_fields(m,l)

     &        ' in list l=', l, ', stat_fname: ', stat_fname(l)

     &    'when trying to add stat_list l=', l

     &    ' Frq=', stat_freq(l), ', stat_fname: ', stat_fname(l)

         DO l=1,nlevels(n),10

          m = MIN(nlevels(n),l+9)

          WRITE(msgBuf,'(A,1P10E10.3)')' interp:  ', (levs(k,n),k=l,m)

         DO l=1,nlevels(n),20

          m = MIN(nlevels(n),l+19)

          WRITE(msgBuf,'(A,20F5.0)') suffix, (levs(k,n),k=l,m)

          WRITE(msgBuf,'(21A)') ' Fields:   ',(' ',flds(l,n),l=nf,m)

        l = 10

          l = l+3

          IF (l.LE.MAX_LEN_MBUF) WRITE(msgBuf(l-2:l),'(I3)') j

     &                 (' ',diagSt_Flds(l,n),l=nf,m)

      INTEGER j,k,l,n,m,nf

           DO l=1,diagSt_nbLists

            DO k=1,diagSt_nbActv(l)

             IF ( iSdiag(k,l).EQ.-iSp ) THEN

               region2fill(j) = MAX(region2fill(j),diagSt_region(j,l))

      INTEGER m, n, j, k, l, bi, bj

           DO l=1,diagSt_nbLists

            DO k=1,diagSt_nbActv(l)

             IF ( iSdiag(k,l).EQ.-iSp ) THEN

               region2fill(j) = MAX(region2fill(j),diagSt_region(j,l))

      INTEGER m, n, j, k, l, bi, bj

        DO l=1,LEN(ccLine)

         ccLine(l:l) = '-'

         l = ILNBLNK(tdiag(n))

         IF ( hdiag(n).NE.0 .AND. l.GE.1 ) THEN

     &                 udiag(n), '|', tdiag(n)(1:l)

         ELSEIF (l.GE.1) THEN

     &                 udiag(n), '|', tdiag(n)(1:l)

        DO l = 1,n-1

         IF ( cdiag(l).EQ.cdiag(n) ) THEN

     &       '1rst (l=', l, ' ), title= ',tdiag(l)

      INTEGER l, n, ncount

        DO l=1,listId

         IF (flag .AND. freq(l) .EQ.freq(listId)

     &            .AND. phase(l).EQ.phase(listId)

     &            .AND. averageFreq(l) .EQ.averageFreq(listId)

     &            .AND. averagePhase(l).EQ.averagePhase(listId)

     &            .AND. averageCycle(l).EQ.averageCycle(listId) ) THEN

          DO k=1,MIN(nActive(l),numperlist)

           IF (flag .AND. jdiag(k,l).GT.0) THEN

            IF (cdiag(mate).EQ.cdiag(jdiag(k,l)) ) THEN

             mdiag(mId,listId) = ABS(idiag(k,l))

      INTEGER k, l

      DO l=1,listId

       IF (flag .AND. freq(l) .EQ. freq(listId)

     &          .AND. phase(l).EQ.phase(listId)

     &          .AND. averageFreq(l) .EQ.averageFreq(listId)

     &          .AND. averagePhase(l).EQ.averagePhase(listId)

     &          .AND. averageCycle(l).EQ.averageCycle(listId) ) THEN

        DO k=1,MIN(nActive(l),numperlist)

         IF (flag .AND. jdiag(k,l).GT.0) THEN

          IF ( cdiag(ndId).EQ.cdiag(jdiag(k,l)) ) THEN

           idiag(mId,listId) = -ABS(idiag(k,l))

           mdiag(mId,listId) = mdiag(k,l)

             DO l=1,averageCycle(ld)

               ii = ipt+l*kdiag(jdiag(md,ld))

               IF (l.EQ.averageCycle(ld)) THEN

      INTEGER j, k, k1, k2, l

        DO l=1,LEN(ccLine)

         ccLine(l:l) = '-'

      DO l=1,diagSt_nbLists

        IF ( diagSt_region(j,l).NE.0 ) THEN

     &       ') in list l=', l

      INTEGER j, k, l

           DO l=1,ld-1

            DO k=1,diagSt_nbActv(l)

              IF ( diagSt_Flds(k,l).EQ.cdiag(-mate) ) inList=.TRUE.

     &                 .OR. (diagSt_region(j,l).LT.diagSt_region(j,ld))

        DO l=1,listId

         IF (flag .AND. diagSt_freq(l) .EQ. diagSt_freq(listId)

     &            .AND. diagSt_phase(l).EQ.diagSt_phase(listId) ) THEN

          DO k=1,MIN(diagSt_nbActv(l),numperlist)

           IF (flag .AND. jSdiag(k,l).GT.0) THEN

            IF (cdiag(mate).EQ.cdiag(jSdiag(k,l)) ) THEN

             mSdiag(mId,listId) = ABS(iSdiag(k,l))

      INTEGER k, l

      DO l=1,listId

       IF (flag .AND. diagSt_freq(l) .EQ. diagSt_freq(listId)

     &          .AND. diagSt_phase(l).EQ.diagSt_phase(listId) ) THEN

        DO k=1,MIN(diagSt_nbActv(l),numperlist)

         IF (flag .AND. jSdiag(k,l).GT.0) THEN

          IF (cdiag(ndId).EQ.cdiag(jSdiag(k,l)) ) THEN

           iSdiag(mId,listId) = -ABS(iSdiag(k,l))

           mSdiag(mId,listId) = mSdiag(k,l)

           do l=0,1

            px_ind(l+1) = x_in(w_ind(i,j)+l)

            py_ind(l+1) = y_in(s_ind(i,j)+l)

           do l=-1,2

            px_ind(l+2) = x_in(w_ind(i,j)+l)

            py_ind(l+2) = y_in(s_ind(i,j)+l)

      integer  i, j, k, l, js, bi, bj, sp, interp_unit