/[MITgcm]/MITgcm/pkg/fizhi/fizhi_turb.F
ViewVC logotype

Diff of /MITgcm/pkg/fizhi/fizhi_turb.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.16 by molod, Wed Jul 28 01:25:07 2004 UTC revision 1.30 by molod, Tue Feb 22 22:51:42 2005 UTC
# Line 98  c--------------------------------------- Line 98  c---------------------------------------
98    
99  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
100  #include "SIZE.h"  #include "SIZE.h"
101  #include "diagnostics_SIZE.h"  #include "DIAGNOSTICS_SIZE.h"
102  #include "diagnostics.h"  #include "DIAGNOSTICS.h"
103  #endif  #endif
104    
105        integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb        integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb
# Line 150  C Set fmu and ed to zero for no backgrou Line 150  C Set fmu and ed to zero for no backgrou
150        _RL  fcctmp(im,jm,nlay)        _RL  fcctmp(im,jm,nlay)
151        _RL tmpdiag(im,jm)        _RL tmpdiag(im,jm)
152        _RL   thtgz(im*jm)        _RL   thtgz(im*jm)
153          logical  diagnostics_is_on
154          external diagnostics_is_on
155    
156        integer nland        integer nland
157        _RL alwcoeff(nchp),blwcoeff(nchp)        _RL alwcoeff(nchp),blwcoeff(nchp)
# Line 231  C Set fmu and ed to zero for no backgrou Line 233  C Set fmu and ed to zero for no backgrou
233        _RL qliqmsc(nchp,nlay),fccmsc(nchp,nlay)        _RL qliqmsc(nchp,nlay),fccmsc(nchp,nlay)
234    
235        integer ndlsm        integer ndlsm
236        parameter ( ndlsm = 40)        parameter ( ndlsm = 1)
237        _RL qdiaglsm(nchp,ndlsm)        _RL qdiaglsm(nchp,ndlsm)
238    
       integer n,nsecf,nmonf,ndayf  
       nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)  
       nmonf(n) = mod(n,10000)/100  
       ndayf(n) = mod(n,100)  
   
239        _RL pi,secday,sdayopi2,rgas,akap,cp,alhl        _RL pi,secday,sdayopi2,rgas,akap,cp,alhl
240        _RL faceps,grav,caltoj,virtcon,getcon        _RL faceps,grav,caltoj,virtcon,getcon
241        _RL heatw,undef,timstp,delttrb,dttrb,ra        _RL heatw,undef,timstp,delttrb,dttrb,ra
# Line 247  C Set fmu and ed to zero for no backgrou Line 244  C Set fmu and ed to zero for no backgrou
244        integer nocean, nice        integer nocean, nice
245        integer ndmoist,time_left,ndum        integer ndmoist,time_left,ndum
246        integer ntracedim        integer ntracedim
247        _RL    dtfac,timstp2,sum        _RL    dtfac,timstp2,sum0
248  C logical begin flag - set to true to indicate a cold start  C logical begin flag - set to true to indicate a cold start
249        logical qbeg            logical qbeg    
250    
251          integer n,nsecf,nmonf,ndayf
252          nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)
253          nmonf(n) = mod(n,10000)/100
254          ndayf(n) = mod(n,100)
255    
256  #ifdef CRAY  #ifdef CRAY
257  #ifdef f77  #ifdef f77
258  cfpp$ expand (qsat)  cfpp$ expand (qsat)
# Line 259  cfpp$ expand (qsat) Line 261  cfpp$ expand (qsat)
261    
262  c   compute variables that do not change  c   compute variables that do not change
263  c  c
264    
265         pi = 4.*atan(1.)         pi = 4.*atan(1.)
266         secday   = getcon('SDAY')         secday   = getcon('SDAY')
267         sdayopi2 = getcon('SDAY') / (pi*2.)         sdayopi2 = getcon('SDAY') / (pi*2.)
# Line 303  c ************************************** Line 306  c **************************************
306                
307        qbeg = .false.        qbeg = .false.
308    
309        sum = 0.0        sum0 = 0.0
310        do L=1,nlay        do L=1,nlay
311        do n=1,nchptot        do n=1,nchptot
312        sum = sum + tke(n,L)        sum0 = sum0 + tke(n,L)
313        enddo        enddo
314        enddo        enddo
315    
316  #ifdef ALLOW_USE_MPI  #ifdef ALLOW_USE_MPI
317        call mpi_allreduce(sum,const,1,mpi_double_precision,mpi_sum,        call mpi_allreduce(sum0,const,1,mpi_double_precision,mpi_sum,
318       .                                                mpi_comm_world,n)       .                                                mpi_comm_world,n)
319  #else  #else
320        const = sum        const = sum0
321  #endif  #endif
322    
323        if( const.eq.0.0 ) then        if( const.eq.0.0 ) then
324        qbeg = .true.        qbeg = .true.
325            if( myid.eq.0 ) then            if( myid.eq.1 .and. bi.eq.1 ) then
326            print *            print *
327            print *, 'Warning!'            print *, 'Warning!'
328            print *, 'Turbulent Kinetic Energy has not been initialized.'            print *, 'Turbulent Kinetic Energy has not been initialized.'
# Line 328  c ************************************** Line 331  c **************************************
331            endif            endif
332        endif        endif
333    
334          call diagnostics_fill(tgz,'TGROUND ',0,1,3,bi,bj,myid)
335    c     if(itground.gt.0) then
336    c     do j =1,jm
337    c     do i =1,im
338    c     qdiag(i,j,itground,bi,bj) = qdiag(i,j,itground,bi,bj) + tgz(i,j)
339    c     enddo
340    c     enddo
341    c     endif
342    
343  c **********************************************************************  c **********************************************************************
344  c                            Initialization  c                            Initialization
345  c **********************************************************************  c **********************************************************************
# Line 335  c ************************************** Line 347  c **************************************
347  c Initialize diagnostic for ground temperature change  c Initialize diagnostic for ground temperature change
348  c ---------------------------------------------------  c ---------------------------------------------------
349        if(idtg.gt.0) then        if(idtg.gt.0) then
350        do i =1,ijall        do j =1,jm
351        qdiag(i,1,idtg,bi,bj) = qdiag(i,1,idtg,bi,bj) - tgz(i,1)        do i =1,im
352          qdiag(i,j,idtg,bi,bj) = qdiag(i,j,idtg,bi,bj) - tgz(i,j)
353          enddo
354        enddo        enddo
355        endif        endif
356    
# Line 348  c ************************************** Line 362  c **************************************
362    
363        numstrips   = ((nchptot-1)/istrip) + 1        numstrips   = ((nchptot-1)/istrip) + 1
364    
       print *,' In turb - about to call grd2msc for ',nchptot,' tiles '  
       print *,' im ',im,' jm ',jm,' nchp ',nchp,' ijall ',ijall  
       print *,' igrd ',(igrd(i),i=1,nchptot)  
       print *,' pz ',(pz(i,1),i=1,ijall)  
       stop  
   
365        call grd2msc(pz(1,1),im,jm,igrd,pmsc,nchp,nchptot)        call grd2msc(pz(1,1),im,jm,igrd,pmsc,nchp,nchptot)
366    
367        call grd2msc(tgz,im,jm,igrd,tground,nchp,nchptot)        call grd2msc(tgz,im,jm,igrd,tground,nchp,nchptot)
# Line 461  c ----------------------------------- Line 469  c -----------------------------------
469        tmpdiag(i,1) = 0.0        tmpdiag(i,1) = 0.0
470        enddo        enddo
471        call msc2grd (igrd,chfr,qice,nchp,nchptot,fracland,tmpdiag,im,jm)        call msc2grd (igrd,chfr,qice,nchp,nchptot,fracland,tmpdiag,im,jm)
472        do i =1,ijall        do j =1,jm
473        qdiag(i,1,iqice,bi,bj) = qdiag(i,1,iqice,bi,bj) + tmpdiag(i,1)        do i =1,im
474          qdiag(i,j,iqice,bi,bj) = qdiag(i,j,iqice,bi,bj) + tmpdiag(i,j)
475          enddo
476        enddo        enddo
477        nqice = nqice + 1        nqice = nqice + 1
478        endif        endif
479    
       print *,' In turb - did all initialization - ready for regions '  
       stop  
480  c**********************************************************************  c**********************************************************************
481  c                        loop over regions  c                        loop over regions
482  c**********************************************************************  c**********************************************************************
# Line 484  c*************************************** Line 492  c***************************************
492         call strip2tile(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)         call strip2tile(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)
493         call strip2tile(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)         call strip2tile(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)
494         call strip2tile(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)         call strip2tile(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)
495         do nt = 1,ntracers-ptracers  c      do nt = 1,ntracers-ptracers
496         call strip2tile(qz(1,1,1,ptracers+nt),igrd,tracers(1,1,nt),nchp,  c      call strip2tile(qz(1,1,1,ptracers+nt),igrd,tracers(1,1,nt),nchp,
497       1                                             ijall,istrip,nlay,nn)  c    1                                             ijall,istrip,nlay,nn)
498         enddo  c      enddo
499    
500         call stripit  (z0,stz0,nchptot,nchp,istrip,1,nn)         call stripit  (z0,stz0,nchptot,nchp,istrip,1,nn)
501         call stripit  (tground,th(1,nlay+1),nchptot,nchp,istrip,1,nn)         call stripit  (tground,th(1,nlay+1),nchptot,nchp,istrip,1,nn)
# Line 567  c -------------------------------------- Line 575  c --------------------------------------
575        sh(i,nlay+1) = qa(i)        sh(i,nlay+1) = qa(i)
576        enddo        enddo
577    
578        if(iqg.gt.0) then  c     if(iqg.gt.0) then
579        do i=1,istrip  c     do i=1,istrip
580        tmpstrip(i) = sh(i,nlay+1)*1000  c     tmpstrip(i) = sh(i,nlay+1)*1000
581        enddo  c     enddo
582        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
583       1                        qdiag(1,1,iqg,bi,bj),ijall,1,nn,.false.)  c    1                        qdiag(1,1,iqg,bi,bj),ijall,1,nn,.false.)
584        endif  c     endif
585    
586  c value of tracers at the ground  c value of tracers at the ground
587  c ------------------------------  c ------------------------------
588        do nt = 1,ntracers-ptracers  c     do nt = 1,ntracers-ptracers
589         do i = 1,istrip  C      do i = 1,istrip
590          tracers(i,nlay+1,nt) = 0.  C       tracers(i,nlay+1,nt) = 0.
591         enddo  C      enddo
592        enddo  C     enddo
593    
594  c compute virtual potential temperatures  c compute virtual potential temperatures
595  c --------------------------------------  c --------------------------------------
# Line 617  c ------------------------------- Line 625  c -------------------------------
625    
626  c increment diagnostic arrays for quantities calculated before trbfl  c increment diagnostic arrays for quantities calculated before trbfl
627  c ------------------------------------------------------------------  c ------------------------------------------------------------------
628        do i =1,istrip  c     do i =1,istrip
629        stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)  c     stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)
630        enddo  c     enddo
631        if(idtsrf.gt.0) then  c     if(idtsrf.gt.0) then
632        call paste2grd(stdiag,igrd,chfrstr,istrip,nchp,  c     call paste2grd(stdiag,igrd,chfrstr,istrip,nchp,
633       1                      qdiag(1,1,idtsrf,bi,bj),ijall,1,nn,.false.)  c    1                      qdiag(1,1,idtsrf,bi,bj),ijall,1,nn,.false.)
634    c     endif
635    
636    
637          if(2.eq.1)then
638          print *,' In turb before trbflx - strip ',nn,' out of ',numstrips
639          print *,' bi = ',bi
640          print *,' ntracers ',ntracers,' ptracers ',ptracers
641          print *,'dttrb,itrtrb,rmu,edle ',dttrb,' ',itrtrb,' ',rmu,' ',edle
642          print *,' nchp ',nchp,' nchptot ',nchptot,' nchplnd ',nchplnd
643          print *,' qbeg, tprof ',qbeg,' ',tprof
644          print *,'istrip,nlay,nymd,nhms ',istrip,' ',nlay,' ',nymd,' ',nhms
645          print *,' grav,cp,rgas,faceps,virtcon,undef ',
646         .     grav,' ',cp,' ',rgas,' ',faceps,' ',virtcon,' ',undef
647          print *,' field: th ',th
648    c     print *,' field: thv ',thv
649    c     print *,' field: sh ',sh
650    c     print *,' field: u ',u
651    c     print *,' field: v ',v
652          print *,' field: p ',p
653    c     print *,' field: pe ',pe
654    c     print *,' field: pk ',pk
655    c     print *,' field: pke ',pke
656    c     print *,' field: dpstr ',dpstr
657    c     print *,' field: stwatr ',stwatr
658    c     print *,' field: stz0 ',stz0
659        endif        endif
660    
661  c call trbflx  c call trbflx
# Line 636  c ----------- Line 669  c -----------
669       6 hsturb,dhsdqa,dhsdtc,transth,transsh,       6 hsturb,dhsdqa,dhsdtc,transth,transsh,
670       7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)       7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
671    
672    
673          if(2.eq.1)then
674          print *,' In turbio, Just after trbflx for strip ',nn,' bi = ',bi
675          print *,' field: stuflux ',stuflux
676          print *,' field: stvflux ',stvflux
677          print *,' field: dshdthg ',dshdthg
678          print *,' field: dshdshg ',dshdshg
679          print *,' field: dthdthg ',dthdthg
680          print *,' field: dthdshg ',dthdshg
681          print *,' field: scu ',scu
682          print *,' field: eturb ',eturb
683          print *,' field: dedqa ',dedqa
684          print *,' field: dedtc ',dedtc
685          print *,' field: hsturb ',hsturb
686          print *,' field: dhsdqa ',dhsdqa
687          print *,' field: dhsdtc ',dhsdtc
688          print *,' field: transth ',transth
689          print *,' field: transsh ',transsh
690          endif
691    
692        call pastit (qq,tke,istrip,nchp,nchptot,nlay,nn)        call pastit (qq,tke,istrip,nchp,nchptot,nlay,nn)
693        call pastit (ctsave,ctmt,istrip,nchp,nchptot,1,nn)        call pastit (ctsave,ctmt,istrip,nchp,nchptot,1,nn)
694        call pastit (xxsave,xxmt,istrip,nchp,nchptot,1,nn)        call pastit (xxsave,xxmt,istrip,nchp,nchptot,1,nn)
# Line 667  C*************************************** Line 720  C***************************************
720         hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)         hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)
721         psurf(i) = pe(i,nlay+1)         psurf(i) = pe(i,nlay+1)
722         wspeed(i) = sqrt(u(i,nlay)*u(i,nlay) + v(i,nlay)*v(i,nlay))         wspeed(i) = sqrt(u(i,nlay)*u(i,nlay) + v(i,nlay)*v(i,nlay))
723           if(wspeed(i) .lt. 1.e-20) wspeed(i) = 1.e-20
724  C Note: This LSM precip bug needs to be cleaned up  C Note: This LSM precip bug needs to be cleaned up
725  ccc   newsnow(i) = newsnow(i)*dtfac    ccc   newsnow(i) = newsnow(i)*dtfac  
726  ccc   raincon(i) = raincon(i)*dtfac    ccc   raincon(i) = raincon(i)*dtfac  
# Line 706  ccc   rainls (i) = rainls (i)*dtfac Line 760  ccc   rainls (i) = rainls (i)*dtfac
760  c**********************************************************************  c**********************************************************************
761  c   diagnostics: fill arrays for lsm input fields  c   diagnostics: fill arrays for lsm input fields
762  c**********************************************************************  c**********************************************************************
763        if(isnowfall.gt.0) then  c     if(isnowfall.gt.0) then
764        do i = 1,istrip  c     do i = 1,istrip
765        tmpstrip(i) = newsnow(i)*86400    c     tmpstrip(i) = newsnow(i)*86400  
766        enddo  c     enddo
767        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
768       1                   qdiag(1,1,isnowfall,bi,bj),ijall,1,nn,.false.)  c    1                   qdiag(1,1,isnowfall,bi,bj),ijall,1,nn,.false.)
769        endif  c     endif
770        if(iraincon.gt.0) then        if(diagnostics_is_on('RAINCON ',myid) ) then
771        do i = 1,istrip        do i = 1,istrip
772        tmpstrip(i) = raincon(i)*86400          tmpstrip(i) = raincon(i)*86400  
773        enddo        enddo
774        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
775       1                   qdiag(1,1,iraincon,bi,bj),ijall,1,nn,.false.)       . .false., 'RAINCON ', 1, 1, bi, bj, myid)
776        endif        endif
777        if(irainlsp.gt.0) then        if(diagnostics_is_on('RAINLSP ',myid) ) then
778        do i = 1,istrip        do i = 1,istrip
779        tmpstrip(i) = rainls(i)*86400          tmpstrip(i) = rainls(i)*86400  
780        enddo        enddo
781        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
782       1                  qdiag(1,1,irainlsp,bi,bj),ijall,1,nn,.false.)       . .false., 'RAINLSP ', 1, 1, bi, bj, myid)
       endif  
       if(igreen.gt.0) then  
       call paste2grd(grnstrip,igrd,chfrstr,istrip,nchp,  
      1                    qdiag(1,1,igreen,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ilai.gt.0) then  
       call paste2grd(laistrip,igrd,chfrstr,istrip,nchp,  
      1                    qdiag(1,1,ilai,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ipardr.gt.0) then  
       call paste2grd(pardr,igrd,chfrstr,istrip,nchp,  
      1                 qdiag(1,1,ipardr,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ipardf.gt.0) then  
       call paste2grd(pardf,igrd,chfrstr,istrip,nchp,  
      1                 qdiag(1,1,ipardf,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(idlwdtc.gt.0) then  
       call paste2grd(blwrad,igrd,chfrstr,istrip,nchp,  
      1                  qdiag(1,1,idlwdtc,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(idhdtc.gt.0) then  
       call paste2grd(dhsdtc,igrd,chfrstr,istrip,nchp,  
      1                  qdiag(1,1,idhdtc,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(idedtc.gt.0) then  
       call paste2grd(dedtc,igrd,chfrstr,istrip,nchp,  
      1                 qdiag(1,1,idedtc,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(idhdqa.gt.0) then  
       call paste2grd(dhsdqa,igrd,chfrstr,istrip,nchp,  
      1                  qdiag(1,1,idhdqa,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(idedqa.gt.0) then  
       call paste2grd(dedqa,igrd,chfrstr,istrip,nchp,  
      1                 qdiag(1,1,idedqa,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ilwgdown.gt.0) then  
       call paste2grd(hlwdwn,igrd,chfrstr,istrip,nchp,  
      1                  qdiag(1,1,ilwgdown,bi,bj),ijall,1,nn,.false.)  
783        endif        endif
784    c     if(igreen.gt.0) then
785    c     call paste2grd(grnstrip,igrd,chfrstr,istrip,nchp,
786    c    1                    qdiag(1,1,igreen,bi,bj),ijall,1,nn,.false.)
787    c     endif
788    c     if(ilai.gt.0) then
789    c     call paste2grd(laistrip,igrd,chfrstr,istrip,nchp,
790    c    1                    qdiag(1,1,ilai,bi,bj),ijall,1,nn,.false.)
791    c     endif
792          if(diagnostics_is_on('PARDR   ',myid) ) then
793          call diag_vegtile_fill (pardr,igrd,chfrstr,istrip,nchp,nn,
794         . .false., 'PARDR   ', 1, 1, bi, bj, myid)
795          endif
796    c     if(ipardf.gt.0) then
797    c     call paste2grd(pardf,igrd,chfrstr,istrip,nchp,
798    c    1                 qdiag(1,1,ipardf,bi,bj),ijall,1,nn,.false.)
799    c     endif
800    c     if(idlwdtc.gt.0) then
801    c     call paste2grd(blwrad,igrd,chfrstr,istrip,nchp,
802    c    1                  qdiag(1,1,idlwdtc,bi,bj),ijall,1,nn,.false.)
803    c     endif
804    c     if(idhdtc.gt.0) then
805    c     call paste2grd(dhsdtc,igrd,chfrstr,istrip,nchp,
806    c    1                  qdiag(1,1,idhdtc,bi,bj),ijall,1,nn,.false.)
807    c     endif
808    c     if(idedtc.gt.0) then
809    c     call paste2grd(dedtc,igrd,chfrstr,istrip,nchp,
810    c    1                 qdiag(1,1,idedtc,bi,bj),ijall,1,nn,.false.)
811    c     endif
812    c     if(idhdqa.gt.0) then
813    c     call paste2grd(dhsdqa,igrd,chfrstr,istrip,nchp,
814    c    1                  qdiag(1,1,idhdqa,bi,bj),ijall,1,nn,.false.)
815    c     endif
816    c     if(idedqa.gt.0) then
817    c     call paste2grd(dedqa,igrd,chfrstr,istrip,nchp,
818    c    1                 qdiag(1,1,idedqa,bi,bj),ijall,1,nn,.false.)
819    c     endif
820    c     if(ilwgdown.gt.0) then
821    c     call paste2grd(hlwdwn,igrd,chfrstr,istrip,nchp,
822    c    1                  qdiag(1,1,ilwgdown,bi,bj),ijall,1,nn,.false.)
823    c     endif
824  c**********************************************************************  c**********************************************************************
825    
826        if(nland.gt.0)then        if(nland.gt.0)then
827    
828          if(2.eq.1)then
829          print *,' In turbio, Just before tile for strip ',nn,' bi = ',bi
830          print *,' calling tile for ',nland,' land points '
831          print *,' field: types ',types
832          print *,' field: chfrstr ',chfrstr
833    c     print *,' field: rainls ',rainls
834    c     print *,' field: newsnow ',newsnow
835    c     print *,' field: wspeed ',wspeed
836          print *,' field: eturb ',eturb
837          print *,' field: dedqa ',dedqa
838          print *,' field: dedtc ',dedtc
839          print *,' field: hsturb ',hsturb
840          print *,' field: dhsdqa ',dhsdqa
841          print *,' field: dhsdtc ',dhsdtc
842    c     print *,' field: tmpnlay ',tmpnlay
843    c     print *,' field: sh(nlay) ',(sh(i,nlay),i=1,istrip)
844    c     print *,' field: cd ',cd
845    c     print *,' field: cosz ',cosz
846    c     print *,' field: pardr ',pardr
847    c     print *,' field: pardf ',pardf
848          print *,' field: swnet ',swnet
849          print *,' field: hlwdwn ',hlwdwn
850    c     print *,' field: psurf ',psurf
851    c     print *,' field: laistrip ',laistrip
852    c     print *,' field: grnstrip ',grnstrip
853    c     print *,' field: z2str ',z2str
854    c     print *,' field: scatstr ',scatstr
855    c     print *,' field: z2str ',z2str
856    c     print *,' field: rs1str ',rs1str
857    c     print *,' field: rs1str ',rs2str
858    c     print *,' field: rdcstr ',rdcstr
859    c     print *,' field: u2fstr ',u2fstr
860    c     print *,' field: dqsdtstr ',dqsdtstr
861          print *,' field: alwrad ',alwrad
862          print *,' field: blwrad ',blwrad
863          print *,' field: tc ',tc
864          print *,' field: td ',td
865    c     print *,' field: swet1 ',swet1
866    c     print *,' field: swet2 ',swet2
867    c     print *,' field: swet3 ',swet3
868    c     print *,' field: capacity ',capacity
869    c     print *,' field: snowdepth ',snowdepth
870          endif
871    
872         call tile (         call tile (
873       I   nland, timstp, types, rainls, raincon,  newsnow,  wspeed,       I   nland, timstp, types, rainls, raincon,  newsnow,  wspeed,
874       I   eturb,  dedqa,  dedtc,  hsturb, dhsdqa, dhsdtc,       I   eturb,  dedqa,  dedtc,  hsturb, dhsdqa, dhsdtc,
# Line 785  c*************************************** Line 884  c***************************************
884       O   strdg5, strdg6, strdg7, strdg8, strdg9)       O   strdg5, strdg6, strdg7, strdg8, strdg9)
885        endif        endif
886    
887          if(2.eq.1)then
888          print *,' In turbio, Just after tile for strip ',nn
889          print *,' calling tile for ',nland,' land points '
890          print *,' field: tc ',tc
891    c     print *,' field: td ',td
892          print *,' field: strdg1 ',strdg1
893          print *,' field: strdg2 ',strdg2
894          print *,' field: strdg3 ',strdg3
895          print *,' field: strdg4 ',strdg4
896          print *,' field: strdg5 ',strdg5
897          print *,' field: strdg6 ',strdg6
898          print *,' field: strdg7 ',strdg7
899          print *,' field: strdg8 ',strdg8
900          print *,' field: strdg9 ',strdg9
901    c     print *,' field: swet1 ',swet1
902    c     print *,' field: swet2 ',swet2
903    c     print *,' field: swet3 ',swet3
904    c     print *,' field: capacity ',capacity
905    c     print *,' field: snowdepth ',snowdepth
906          endif
907        if( nice.gt.0 ) then        if( nice.gt.0 ) then
908           print *,' Calling sea ice routine - SHOULD NOT BE HERE!'
909         call seaice ( nocean, timstp,     hice,         call seaice ( nocean, timstp,     hice,
910       .               eturb(nland+1),    dedtc(nland+1),       .               eturb(nland+1),    dedtc(nland+1),
911       .              hsturb(nland+1),   dhsdtc(nland+1),       .              hsturb(nland+1),   dhsdtc(nland+1),
# Line 799  c*************************************** Line 919  c***************************************
919  c   diagnostics: fill arrays for lsm output fields  c   diagnostics: fill arrays for lsm output fields
920  c***********************************************************************  c***********************************************************************
921    
922        if(irunoff.gt.0) then  c     if(irunoff.gt.0) then
923        call paste2grd(runoff,igrd,chfrstr,istrip,nchp,  c     call paste2grd(runoff,igrd,chfrstr,istrip,nchp,
924       1                  qdiag(1,1,irunoff,bi,bj),ijall,1,nn,.false.)  c    1                  qdiag(1,1,irunoff,bi,bj),ijall,1,nn,.false.)
925        endif  c     endif
926        if(ifwsoil.gt.0) then  c     if(ifwsoil.gt.0) then
927        call paste2grd(fwsoil,igrd,chfrstr,istrip,nchp,  c     call paste2grd(fwsoil,igrd,chfrstr,istrip,nchp,
928       1                  qdiag(1,1,ifwsoil,bi,bj),ijall,1,nn,.false.)  c    1                  qdiag(1,1,ifwsoil,bi,bj),ijall,1,nn,.false.)
929        endif  c     endif
930        if(igdrain.gt.0) then  c     if(igdrain.gt.0) then
931        call paste2grd(gdrain,igrd,chfrstr,istrip,nchp,  c     call paste2grd(gdrain,igrd,chfrstr,istrip,nchp,
932       1                  qdiag(1,1,igdrain,bi,bj),ijall,1,nn,.false.)  c    1                  qdiag(1,1,igdrain,bi,bj),ijall,1,nn,.false.)
933        endif  c     endif
934        if(isnowmelt.gt.0) then  c     if(isnowmelt.gt.0) then
935        call paste2grd(smelt,igrd,chfrstr,istrip,nchp,  c     call paste2grd(smelt,igrd,chfrstr,istrip,nchp,
936       1                 qdiag(1,1,isnowmelt,bi,bj),ijall,1,nn,.false.)  c    1                 qdiag(1,1,isnowmelt,bi,bj),ijall,1,nn,.false.)
937        endif  c     endif
938        if(ieveg.gt.0) then  c     if(ieveg.gt.0) then
939        call paste2grd(eveg,igrd,chfrstr,istrip,nchp,  c     call paste2grd(eveg,igrd,chfrstr,istrip,nchp,
940       1                    qdiag(1,1,ieveg,bi,bj),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ieveg,bi,bj),ijall,1,nn,.false.)
941        endif  c     endif
942        if(iesnow.gt.0) then  c     if(iesnow.gt.0) then
943        call paste2grd(esno,igrd,chfrstr,istrip,nchp,  c     call paste2grd(esno,igrd,chfrstr,istrip,nchp,
944       1                    qdiag(1,1,iesnow,bi,bj),ijall,1,nn,.false.)  c    1                    qdiag(1,1,iesnow,bi,bj),ijall,1,nn,.false.)
945        endif  c     endif
946        if(iesoil.gt.0) then  c     if(iesoil.gt.0) then
947        call paste2grd(esoi,igrd,chfrstr,istrip,nchp,  c     call paste2grd(esoi,igrd,chfrstr,istrip,nchp,
948       1                    qdiag(1,1,iesoil,bi,bj),ijall,1,nn,.false.)  c    1                    qdiag(1,1,iesoil,bi,bj),ijall,1,nn,.false.)
949        endif  c     endif
950        if(ieresv.gt.0) then  c     if(ieresv.gt.0) then
951        call paste2grd(eint,igrd,chfrstr,istrip,nchp,  c     call paste2grd(eint,igrd,chfrstr,istrip,nchp,
952       1                    qdiag(1,1,ieresv,bi,bj),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ieresv,bi,bj),ijall,1,nn,.false.)
953        endif  c     endif
954        if(ievpot.gt.0) then  c     if(ievpot.gt.0) then
955        call paste2grd(evpot,igrd,chfrstr,istrip,nchp,  c     call paste2grd(evpot,igrd,chfrstr,istrip,nchp,
956       1                     qdiag(1,1,ievpot,bi,bj),ijall,1,nn,.false.)  c    1                     qdiag(1,1,ievpot,bi,bj),ijall,1,nn,.false.)
957        endif  c     endif
958        if(idtc.gt.0) then  c     if(idtc.gt.0) then
959        call paste2grd(strdg1,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg1,igrd,chfrstr,istrip,nchp,
960       1                      qdiag(1,1,idtc,bi,bj),ijall,1,nn,.false.)  c    1                      qdiag(1,1,idtc,bi,bj),ijall,1,nn,.false.)
961        endif  c     endif
962        if(idqc.gt.0) then  c     if(idqc.gt.0) then
963        call paste2grd(strdg2,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg2,igrd,chfrstr,istrip,nchp,
964       1                      qdiag(1,1,idqc,bi,bj),ijall,1,nn,.false.)  c    1                      qdiag(1,1,idqc,bi,bj),ijall,1,nn,.false.)
965        endif  c     endif
966        if(itcdtc.gt.0) then  c     if(itcdtc.gt.0) then
967        call paste2grd(strdg3,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg3,igrd,chfrstr,istrip,nchp,
968       1                      qdiag(1,1,itcdtc,bi,bj),ijall,1,nn,.false.)  c    1                      qdiag(1,1,itcdtc,bi,bj),ijall,1,nn,.false.)
969        endif  c     endif
970        if(iraddtc.gt.0) then  c     if(iraddtc.gt.0) then
971        call paste2grd(strdg4,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg4,igrd,chfrstr,istrip,nchp,
972       1                      qdiag(1,1,iraddtc,bi,bj),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iraddtc,bi,bj),ijall,1,nn,.false.)
973        endif  c     endif
974        if(isensdtc.gt.0) then  c     if(isensdtc.gt.0) then
975        call paste2grd(strdg5,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg5,igrd,chfrstr,istrip,nchp,
976       1                     qdiag(1,1,isensdtc,bi,bj),ijall,1,nn,.false.)  c    1                     qdiag(1,1,isensdtc,bi,bj),ijall,1,nn,.false.)
977        endif  c     endif
978        if(ilatdtc.gt.0) then  c     if(ilatdtc.gt.0) then
979        call paste2grd(strdg6,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg6,igrd,chfrstr,istrip,nchp,
980       1                      qdiag(1,1,ilatdtc,bi,bj),ijall,1,nn,.false.)  c    1                      qdiag(1,1,ilatdtc,bi,bj),ijall,1,nn,.false.)
981        endif  c     endif
982        if(itddtc.gt.0) then  c     if(itddtc.gt.0) then
983        call paste2grd(strdg7,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg7,igrd,chfrstr,istrip,nchp,
984       1                      qdiag(1,1,itddtc,bi,bj),ijall,1,nn,.false.)  c    1                      qdiag(1,1,itddtc,bi,bj),ijall,1,nn,.false.)
985        endif  c     endif
986        if(iqcdtc.gt.0) then  c     if(iqcdtc.gt.0) then
987        call paste2grd(strdg8,igrd,chfrstr,istrip,nchp,  c     call paste2grd(strdg8,igrd,chfrstr,istrip,nchp,
988       1                      qdiag(1,1,iqcdtc,bi,bj),ijall,1,nn,.false.)  c    1                      qdiag(1,1,iqcdtc,bi,bj),ijall,1,nn,.false.)
989        endif  c     endif
990  c***********************************************************************  c***********************************************************************
991    
992        if( ndlsm.gt.1 ) then        if( ndlsm.gt.1 ) then
# Line 948  c*************************************** Line 1068  c***************************************
1068        enddo        enddo
1069        enddo        enddo
1070    
1071          if(2.eq.1)then
1072          print *,' In turbio, just after updating th and sh - strip ',nn
1073          print *,' field: th ',th
1074          print *,' field: sh ',sh
1075          endif
1076    
1077  c tendency updates  c tendency updates
1078  c ----------------  c ----------------
1079        do  l=1,nlay        do  l=1,nlay
# Line 970  c ---------------- Line 1096  c ----------------
1096        do i =1,istrip        do i =1,istrip
1097        tends(i) = ( th(i,l)-tmpstrip(i) )        tends(i) = ( th(i,l)-tmpstrip(i) )
1098        enddo        enddo
1099    
1100          if(2.eq.1)then
1101          print *,' In turbio, tendencies for strip ',nn,' level ',l
1102          print *,' field: th ',tends
1103          endif
1104    
1105        call pastit (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)        call pastit (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)
1106    
1107        call strip2tile(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,        call strip2tile(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,
# Line 977  c ---------------- Line 1109  c ----------------
1109        do i =1,istrip        do i =1,istrip
1110        tends(i) = ( sh(i,l)-tmpstrip(i) )        tends(i) = ( sh(i,l)-tmpstrip(i) )
1111        enddo        enddo
1112    
1113          if(2.eq.1)then
1114          print *,' In turbio, tendencies for strip ',nn,' level ',l
1115          print *,' field: sh ',tends
1116          endif
1117    
1118        call pastit (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)        call pastit (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)
1119    
1120        do nt = 1,ntracers-ptracers  c     do nt = 1,ntracers-ptracers
1121         call strip2tile(qz(1,1,L,ptracers+nt),igrd,tmpstrip,nchp,  c      call strip2tile(qz(1,1,L,ptracers+nt),igrd,tmpstrip,nchp,
1122       1                                             ijall,istrip,1,nn)  c    1                                             ijall,istrip,1,nn)
1123         do i =1,istrip  c      do i =1,istrip
1124          tends(i) = ( tracers(i,L,nt)-tmpstrip(i) )  c       tends(i) = ( tracers(i,L,nt)-tmpstrip(i) )
1125         enddo  c      enddo
1126         call pastit(tends,dqmsc(1,L,ptracers+nt),istrip,nchp,  c      call pastit(tends,dqmsc(1,L,ptracers+nt),istrip,nchp,
1127       .                                     nchptot,1,nn)  c    .                                     nchptot,1,nn)
1128        enddo  c     enddo
1129    
1130        enddo        enddo
1131    
# Line 999  c note: the order, logic, and scaling of Line 1137  c note: the order, logic, and scaling of
1137  c       diagnostics is critical!  c       diagnostics is critical!
1138  c ------------------------------  c ------------------------------
1139    
1140        if(ievap.gt.0) then  c     if(ievap.gt.0) then
1141        do i=1,istrip  c     do i=1,istrip
1142        tmpstrip(i) = stqflux(i,nlay) * 86400  c     tmpstrip(i) = stqflux(i,nlay) * 86400
1143        enddo  c     enddo
1144        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1145       1                    qdiag(1,1,ievap,bi,bj),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ievap,bi,bj),ijall,1,nn,.false.)
1146        endif  c     endif
1147        if(ieflux.gt.0) then  c     if(ieflux.gt.0) then
1148        do i=1,istrip  c     do i=1,istrip
1149        tmpstrip(i) = stqflux(i,nlay) * alhl  c     tmpstrip(i) = stqflux(i,nlay) * alhl
1150        enddo  c     enddo
1151        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1152       1                    qdiag(1,1,ieflux,bi,bj),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ieflux,bi,bj),ijall,1,nn,.false.)
1153        endif  c     endif
1154        if(ihflux.gt.0) then  c     if(ihflux.gt.0) then
1155        do i=1,istrip  c     do i=1,istrip
1156        tmpstrip(i) = sttflux(i,nlay) * cp  c     tmpstrip(i) = sttflux(i,nlay) * cp
1157        enddo  c     enddo
1158        call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1159       1                    qdiag(1,1,ihflux,bi,bj),ijall,1,nn,.false.)  c    1                    qdiag(1,1,ihflux,bi,bj),ijall,1,nn,.false.)
1160        endif  c     endif
1161        if(ituflux.gt.0) then        if(diagnostics_is_on('TUFLUX  ',myid) ) then
1162        call paste2grd(stuflux,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stuflux,igrd,chfrstr,istrip,nchp,nn,
1163       1                   qdiag(1,1,ituflux,bi,bj),ijall,nlay,nn,.false.)       . .false., 'TUFLUX  ', 0, nlay, bi, bj, myid)
1164        endif        endif
1165        if(itvflux.gt.0) then        if(diagnostics_is_on('TVFLUX  ',myid) ) then
1166        call paste2grd(stvflux,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stvflux,igrd,chfrstr,istrip,nchp,nn,
1167       1                   qdiag(1,1,itvflux,bi,bj),ijall,nlay,nn,.false.)       . .false., 'TVFLUX  ', 0, nlay, bi, bj, myid)
1168        endif        endif
1169        if(ittflux.gt.0) then  c     if(ituflux.gt.0) then
1170    c     call paste2grd(stuflux,igrd,chfrstr,istrip,nchp,
1171    c    1                   qdiag(1,1,ituflux,bi,bj),ijall,nlay,nn,.false.)
1172    c     endif
1173          if(diagnostics_is_on('TTFLUX  ',myid) ) then
1174        do l=1,nlay        do l=1,nlay
1175        do i=1,istrip        do i=1,istrip
1176        sttflux(i,l) = sttflux(i,l) * cp        sttflux(i,l) = sttflux(i,l) * cp
1177        enddo        enddo
1178        enddo        enddo
1179        call paste2grd(sttflux,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (sttflux,igrd,chfrstr,istrip,nchp,nn,
1180       1                   qdiag(1,1,ittflux,bi,bj),ijall,nlay,nn,.false.)       . .false., 'TTFLUX  ', 0, nlay, bi, bj, myid)
1181        endif        endif
1182        if(itqflux.gt.0) then        if(diagnostics_is_on('TQFLUX  ',myid) ) then
1183        do l=1,nlay        do l=1,nlay
1184        do i=1,istrip        do i=1,istrip
1185        stqflux(i,l) = stqflux(i,l) * alhl        stqflux(i,l) = stqflux(i,l) * alhl
1186        enddo        enddo
1187        enddo        enddo
1188        call paste2grd(stqflux,igrd,chfrstr,istrip,nchp,        call diag_vegtile_fill (stqflux,igrd,chfrstr,istrip,nchp,nn,
1189       1                   qdiag(1,1,itqflux,bi,bj),ijall,nlay,nn,.false.)       . .false., 'TQFLUX  ', 0, nlay, bi, bj, myid)
       endif  
       if(iri.gt.0) call paste2grd(sri,igrd,chfrstr,istrip,nchp,  
      1                       qdiag(1,1,iri,bi,bj),ijall,nlay,nn,.false.)  
       if(ikh.gt.0) call paste2grd(skh,igrd,chfrstr,istrip,nchp,  
      1                       qdiag(1,1,ikh,bi,bj),ijall,nlay,nn,.false.)  
       if(ikm.gt.0) call paste2grd(skm,igrd,chfrstr,istrip,nchp,  
      1                       qdiag(1,1,ikm,bi,bj),ijall,nlay,nn,.false.)  
       if(ict.gt.0) then  
       call paste2grd(sct,igrd,chfrstr,istrip,nchp,  
      1                   qdiag(1,1,ict,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(icu.gt.0) then  
       call paste2grd(scu,igrd,chfrstr,istrip,nchp,  
      1                   qdiag(1,1,icu,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(iwinds.gt.0) then  
       call paste2grd(swinds,igrd,chfrstr,istrip,nchp,  
      1                      qdiag(1,1,iwinds,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(iuflux.gt.0) then  
       call paste2grd(stuflux(1,nlay),igrd,chfrstr,istrip,nchp,  
      1                       qdiag(1,1,iuflux,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ivflux.gt.0) then  
       call paste2grd(stvflux(1,nlay),igrd,chfrstr,istrip,nchp,  
      1                       qdiag(1,1,ivflux,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(iustar.gt.0) then  
       call paste2grd(sustar,igrd,chfrstr,istrip,nchp,  
      1                      qdiag(1,1,iustar,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(iz0.gt.0) then  
       call paste2grd(sz0,igrd,chfrstr,istrip,nchp,  
      1                   qdiag(1,1,iz0,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ifrqtrb.gt.0) then  
       call paste2grd(frqtrb,igrd,chfrstr,istrip,nchp,  
      1                      qdiag(1,1,ifrqtrb,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(ipbl.gt.0) then  
       call paste2grd(pbldpth,igrd,chfrstr,istrip,nchp,  
      1                       qdiag(1,1,ipbl,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(iu2m.gt.0) then  
       call paste2grd(stu2m,igrd,chfrstr,istrip,nchp,  
      1                     qdiag(1,1,iu2m,bi,bj),ijall,1,nn,.true.)  
       endif  
       if(iv2m.gt.0) then  
       call paste2grd(stv2m,igrd,chfrstr,istrip,nchp,  
      1                     qdiag(1,1,iv2m,bi,bj),ijall,1,nn,.true.)  
       endif  
       if(it2m.gt.0) then  
       call paste2grd(stt2m,igrd,chfrstr,istrip,nchp,  
      1                     qdiag(1,1,it2m,bi,bj),ijall,1,nn,.true.)  
       endif  
       if(iq2m.gt.0) then  
       do i=1,istrip  
          if( stq2m(i).ne.undef ) then  
              tmpstrip(i) = stq2m(i) * 1000  
          else  
              tmpstrip(i) = undef  
          endif  
       enddo  
       call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  
      1                     qdiag(1,1,iq2m,bi,bj),ijall,1,nn,.true.)  
       endif  
       if(iu10m.gt.0) then  
       call paste2grd(stu10m,igrd,chfrstr,istrip,nchp,  
      1                      qdiag(1,1,iu10m,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(iv10m.gt.0) then  
       call paste2grd(stv10m,igrd,chfrstr,istrip,nchp,  
      1                      qdiag(1,1,iv10m,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(it10m.gt.0) then  
       call paste2grd(stt10m,igrd,chfrstr,istrip,nchp,  
      1                      qdiag(1,1,it10m,bi,bj),ijall,1,nn,.false.)  
       endif  
       if(iq10m.gt.0) then  
       do i=1,istrip  
          if( stq10m(i).ne.undef ) then  
              tmpstrip(i) = stq10m(i) * 1000  
          else  
              tmpstrip(i) = undef  
          endif  
       enddo  
       call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,  
      1                      qdiag(1,1,iq10m,bi,bj),ijall,1,nn,.false.)  
1190        endif        endif
1191          if(diagnostics_is_on('RI      ',myid) ) then
1192          call diag_vegtile_fill (sri,igrd,chfrstr,istrip,nchp,nn,
1193         . .false., 'RI      ', 0, nlay, bi, bj, myid)
1194          endif
1195          if(diagnostics_is_on('KH      ',myid) ) then
1196          call diag_vegtile_fill (skh,igrd,chfrstr,istrip,nchp,nn,
1197         . .false., 'KH      ', 0, nlay, bi, bj, myid)
1198          endif
1199          if(diagnostics_is_on('KM      ',myid) ) then
1200          call diag_vegtile_fill (skm,igrd,chfrstr,istrip,nchp,nn,
1201         . .false., 'KM      ', 0, nlay, bi, bj, myid)
1202          endif
1203    c     if(iri.gt.0) call paste2grd(sri,igrd,chfrstr,istrip,nchp,
1204    c    1                       qdiag(1,1,iri,bi,bj),ijall,nlay,nn,.false.)
1205    c     if(ikh.gt.0) call paste2grd(skh,igrd,chfrstr,istrip,nchp,
1206    c    1                       qdiag(1,1,ikh,bi,bj),ijall,nlay,nn,.false.)
1207    c     if(ikm.gt.0) call paste2grd(skm,igrd,chfrstr,istrip,nchp,
1208    c    1                       qdiag(1,1,ikm,bi,bj),ijall,nlay,nn,.false.)
1209    c     if(ict.gt.0) then
1210    c     call paste2grd(sct,igrd,chfrstr,istrip,nchp,
1211    c    1                   qdiag(1,1,ict,bi,bj),ijall,1,nn,.false.)
1212    c     endif
1213    c     if(icu.gt.0) then
1214    c     call paste2grd(scu,igrd,chfrstr,istrip,nchp,
1215    c    1                   qdiag(1,1,icu,bi,bj),ijall,1,nn,.false.)
1216    c     endif
1217    c     if(iwinds.gt.0) then
1218    c     call paste2grd(swinds,igrd,chfrstr,istrip,nchp,
1219    c    1                      qdiag(1,1,iwinds,bi,bj),ijall,1,nn,.false.)
1220    c     endif
1221    c     if(iuflux.gt.0) then
1222    c     call paste2grd(stuflux(1,nlay),igrd,chfrstr,istrip,nchp,
1223    c    1                       qdiag(1,1,iuflux,bi,bj),ijall,1,nn,.false.)
1224    c     endif
1225    c     if(ivflux.gt.0) then
1226    c     call paste2grd(stvflux(1,nlay),igrd,chfrstr,istrip,nchp,
1227    c    1                       qdiag(1,1,ivflux,bi,bj),ijall,1,nn,.false.)
1228    c     endif
1229    c     if(iustar.gt.0) then
1230    c     call paste2grd(sustar,igrd,chfrstr,istrip,nchp,
1231    c    1                      qdiag(1,1,iustar,bi,bj),ijall,1,nn,.false.)
1232    c     endif
1233    c     if(iz0.gt.0) then
1234    c     call paste2grd(sz0,igrd,chfrstr,istrip,nchp,
1235    c    1                   qdiag(1,1,iz0,bi,bj),ijall,1,nn,.false.)
1236    c     endif
1237    c     if(ifrqtrb.gt.0) then
1238    c     call paste2grd(frqtrb,igrd,chfrstr,istrip,nchp,
1239    c    1                      qdiag(1,1,ifrqtrb,bi,bj),ijall,1,nn,.false.)
1240    c     endif
1241    c     if(ipbl.gt.0) then
1242    c     call paste2grd(pbldpth,igrd,chfrstr,istrip,nchp,
1243    c    1                       qdiag(1,1,ipbl,bi,bj),ijall,1,nn,.false.)
1244    c     endif
1245    c     if(iu2m.gt.0) then
1246    c     call paste2grd(stu2m,igrd,chfrstr,istrip,nchp,
1247    c    1                     qdiag(1,1,iu2m,bi,bj),ijall,1,nn,.true.)
1248    c     endif
1249    c     if(iv2m.gt.0) then
1250    c     call paste2grd(stv2m,igrd,chfrstr,istrip,nchp,
1251    c    1                     qdiag(1,1,iv2m,bi,bj),ijall,1,nn,.true.)
1252    c     endif
1253    c     if(it2m.gt.0) then
1254    c     call paste2grd(stt2m,igrd,chfrstr,istrip,nchp,
1255    c    1                     qdiag(1,1,it2m,bi,bj),ijall,1,nn,.true.)
1256    c     endif
1257    c     if(iq2m.gt.0) then
1258    c     do i=1,istrip
1259    c        if( stq2m(i).ne.undef ) then
1260    c            tmpstrip(i) = stq2m(i) * 1000
1261    c        else
1262    c            tmpstrip(i) = undef
1263    c        endif
1264    c     enddo
1265    c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1266    c    1                     qdiag(1,1,iq2m,bi,bj),ijall,1,nn,.true.)
1267    c     endif
1268    c     if(iu10m.gt.0) then
1269    c     call paste2grd(stu10m,igrd,chfrstr,istrip,nchp,
1270    c    1                      qdiag(1,1,iu10m,bi,bj),ijall,1,nn,.false.)
1271    c     endif
1272    c     if(iv10m.gt.0) then
1273    c     call paste2grd(stv10m,igrd,chfrstr,istrip,nchp,
1274    c    1                      qdiag(1,1,iv10m,bi,bj),ijall,1,nn,.false.)
1275    c     endif
1276    c     if(it10m.gt.0) then
1277    c     call paste2grd(stt10m,igrd,chfrstr,istrip,nchp,
1278    c    1                      qdiag(1,1,it10m,bi,bj),ijall,1,nn,.false.)
1279    c     endif
1280    c     if(iq10m.gt.0) then
1281    c     do i=1,istrip
1282    c        if( stq10m(i).ne.undef ) then
1283    c            tmpstrip(i) = stq10m(i) * 1000
1284    c        else
1285    c            tmpstrip(i) = undef
1286    c        endif
1287    c     enddo
1288    c     call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1289    c    1                      qdiag(1,1,iq10m,bi,bj),ijall,1,nn,.false.)
1290    c     endif
1291    
1292  c**********************************************************************  c**********************************************************************
1293  c   more diagnostics: land surface model parameters  c   more diagnostics: land surface model parameters
1294  c**********************************************************************  c**********************************************************************
1295    
1296        if(itdeep.gt.0)call paste2grd(td,igrd,chfrstr,istrip,nchp,  c     if(itdeep.gt.0)call paste2grd(td,igrd,chfrstr,istrip,nchp,
1297       .                     qdiag(1,1,itdeep,bi,bj),ijall,1,nn,.false.)  c    .                     qdiag(1,1,itdeep,bi,bj),ijall,1,nn,.false.)
1298        if(iqcanopy .gt.0)call paste2grd(qa,igrd,chfrstr,istrip,nchp,  c     if(iqcanopy .gt.0)call paste2grd(qa,igrd,chfrstr,istrip,nchp,
1299       .                  qdiag(1,1,iqcanopy,bi,bj) ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,iqcanopy,bi,bj) ,ijall,1,nn,.false.)
1300        if(ismshal  .gt.0)call paste2grd(swet1,igrd,chfrstr,istrip,nchp,  c     if(ismshal  .gt.0)call paste2grd(swet1,igrd,chfrstr,istrip,nchp,
1301       .                  qdiag(1,1,ismshal,bi,bj)  ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,ismshal,bi,bj)  ,ijall,1,nn,.false.)
1302        if(ismroot  .gt.0)call paste2grd(swet2,igrd,chfrstr,istrip,nchp,  c     if(ismroot  .gt.0)call paste2grd(swet2,igrd,chfrstr,istrip,nchp,
1303       .                  qdiag(1,1,ismroot,bi,bj)  ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,ismroot,bi,bj)  ,ijall,1,nn,.false.)
1304        if(ismdeep  .gt.0)call paste2grd(swet3,igrd,chfrstr,istrip,nchp,  c     if(ismdeep  .gt.0)call paste2grd(swet3,igrd,chfrstr,istrip,nchp,
1305       .                  qdiag(1,1,ismdeep,bi,bj)  ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,ismdeep,bi,bj)  ,ijall,1,nn,.false.)
1306        if(icapacity.gt.0)call paste2grd(capacity,igrd,chfrstr,istrip,  c     if(icapacity.gt.0)call paste2grd(capacity,igrd,chfrstr,istrip,
1307       .             nchp,qdiag(1,1,icapacity,bi,bj),ijall,1,nn,.false.)  c    .             nchp,qdiag(1,1,icapacity,bi,bj),ijall,1,nn,.false.)
1308        if(isnow.gt.0)call paste2grd(snowdepth,igrd,chfrstr,istrip,nchp,  c     if(isnow.gt.0)call paste2grd(snowdepth,igrd,chfrstr,istrip,nchp,
1309       .                  qdiag(1,1,isnow,bi,bj)    ,ijall,1,nn,.false.)  c    .                  qdiag(1,1,isnow,bi,bj)    ,ijall,1,nn,.false.)
1310    
1311  c**********************************************************************  c**********************************************************************
1312  c end regions loop  c end regions loop
# Line 1270  c **** increment diagnostic array for gr Line 1424  c **** increment diagnostic array for gr
1424  c ***       ground temp tendency, and ground humidity  c ***       ground temp tendency, and ground humidity
1425  c *********************************************************************  c *********************************************************************
1426    
1427        if(itground.gt.0) then  c     if(itground.gt.0) then
1428        do i =1,ijall  c     do j =1,jm
1429        qdiag(i,1,itground,bi,bj) = qdiag(i,1,itground,bi,bj) + tgz(i,1)  c     do i =1,im
1430        enddo  c     qdiag(i,j,itground,bi,bj) = qdiag(i,j,itground,bi,bj) + tgz(i,j)
1431        endif  c     enddo
1432    c     enddo
1433    c     endif
1434    
1435        if(itcanopy.gt.0) then        if(itcanopy.gt.0) then
1436        do i =1,ijall        do j =1,jm
1437        qdiag(i,1,itcanopy,bi,bj) = qdiag(i,1,itcanopy,bi,bj) + tgz(i,1)        do i =1,im
1438          qdiag(i,j,itcanopy,bi,bj) = qdiag(i,j,itcanopy,bi,bj) + tgz(i,j)
1439          enddo
1440        enddo        enddo
1441        endif        endif
1442    
1443        if(its.gt.0) then        if(its.gt.0) then
1444        do i =1,ijall        do j =1,jm
1445        tmpstrip(i) = tz(i,1,nlay) * pkht(i,1,nlay)        do i =1,im
1446          tmpdiag(i,j) = tz(i,j,nlay) * pkht(i,j,nlay)
1447          enddo
1448          enddo
1449          do j =1,jm
1450          do i =1,im
1451          qdiag(i,j,its,bi,bj) = qdiag(i,j,its,bi,bj) + tmpdiag(i,j)
1452        enddo        enddo
       do i =1,ijall  
       qdiag(i,1,its,bi,bj) = qdiag(i,1,its,bi,bj) + tmpstrip(i)  
1453        enddo        enddo
1454        endif        endif
1455    
1456        if(idtg.gt.0) then        if(idtg.gt.0) then
1457        do i =1,ijall        do j =1,jm
1458        qdiag(i,1,idtg,bi,bj) = qdiag(i,1,idtg,bi,bj) + tgz(i,1)        do i =1,im
1459          qdiag(i,j,idtg,bi,bj) = qdiag(i,j,idtg,bi,bj) + tgz(i,j)
1460          enddo
1461        enddo        enddo
1462        endif        endif
1463    
# Line 1303  c ************************************** Line 1467  c **************************************
1467        do L = 1,nlay        do L = 1,nlay
1468    
1469         if(iturbu.gt.0) then         if(iturbu.gt.0) then
1470         do i =1,ijall         do j =1,jm
1471          qdiag(i,1,iturbu+l-1,bi,bj) =  qdiag(i,1,iturbu+l-1,bi,bj)         do i =1,im
1472       .                      + duturb(i,1,l) * atimstp * secday          qdiag(i,j,iturbu+l-1,bi,bj) =  qdiag(i,j,iturbu+l-1,bi,bj)
1473         .                      + duturb(i,j,l) * atimstp * secday
1474           enddo
1475         enddo         enddo
1476         endif         endif
1477    
1478         if(iturbv.gt.0) then         if(iturbv.gt.0) then
1479         do i =1,ijall         do j =1,jm
1480          qdiag(i,1,iturbv+l-1,bi,bj) =  qdiag(i,1,iturbv+l-1,bi,bj)         do i =1,im
1481       .                      + dvturb(i,1,l) * atimstp * secday          qdiag(i,j,iturbv+l-1,bi,bj) =  qdiag(i,j,iturbv+l-1,bi,bj)
1482         .                      + dvturb(i,j,l) * atimstp * secday
1483           enddo
1484         enddo         enddo
1485         endif         endif
1486    
1487         if(iturbq.gt.0) then         if(iturbq.gt.0) then
1488         do i =1,ijall         do j =1,jm
1489          qdiag(i,1,iturbq+l-1,bi,bj) =  qdiag(i,1,iturbq+l-1,bi,bj)         do i =1,im
1490       .                      + dqturb(i,1,l,1) * atimstp * secday * 1000          qdiag(i,j,iturbq+l-1,bi,bj) =  qdiag(i,j,iturbq+l-1,bi,bj)
1491         .                      + dqturb(i,j,l,1) * atimstp * secday * 1000
1492           enddo
1493         enddo         enddo
1494         endif         endif
1495    
1496         if(iturbt.gt.0) then         if(iturbt.gt.0) then
1497         do i =1,ijall         do j =1,jm
1498          qdiag(i,1,iturbt+l-1,bi,bj) =  qdiag(i,1,iturbt+l-1,bi,bj)         do i =1,im
1499       .                   + dtturb(i,1,l) * pkz(i,1,l)*atimstp*secday          qdiag(i,j,iturbt+l-1,bi,bj) =  qdiag(i,j,iturbt+l-1,bi,bj)
1500         .                   + dtturb(i,j,l) * pkz(i,j,l)*atimstp*secday
1501           enddo
1502         enddo         enddo
1503         endif         endif
1504    
# Line 1368  c ************************************** Line 1540  c **************************************
1540  c ****                bump diagnostic counters                      ***  c ****                bump diagnostic counters                      ***
1541  c *********************************************************************  c *********************************************************************
1542    
1543    #ifdef ALLOW_DIAGNOSTICS
1544          if( (bi.eq.1) .and. (bj.eq.1) ) then
1545        nturbt    = nturbt    + 1        nturbt    = nturbt    + 1
1546        nturbq    = nturbq    + 1        nturbq    = nturbq    + 1
1547        nturbu    = nturbu    + 1        nturbu    = nturbu    + 1
# Line 1446  c ************************************** Line 1620  c **************************************
1620        ntrbqliq  = ntrbqliq  + 1        ntrbqliq  = ntrbqliq  + 1
1621        ntrbfcc   = ntrbfcc   + 1        ntrbfcc   = ntrbfcc   + 1
1622    
1623          nsdiag1 = nsdiag1 + 1
1624          nsdiag2 = nsdiag2 + 1
1625          endif
1626    
1627    #endif
1628    
1629        return        return
1630        end        end
1631        SUBROUTINE TRBFLX (NN,TH,THV,SH,U,V,QQ,PL,PLE,PLK,PLKE,DPSTR,        SUBROUTINE TRBFLX (NN,TH,THV,SH,U,V,QQ,PL,PLE,PLK,PLKE,DPSTR,
# Line 1481  C    Z0            -         SURFACE ROU Line 1661  C    Z0            -         SURFACE ROU
1661  C    tracers       -         array of passive tracers  C    tracers       -         array of passive tracers
1662  C    ntrace        -         number of tracers to be diffused  C    ntrace        -         number of tracers to be diffused
1663  C    ntracedim     -         outer dimension of tracers array  C    ntracedim     -         outer dimension of tracers array
1664  C    DTAU          -         TIME CHANGE PER ITERATION OF TRBLFX  C    DTAU          -         TIME CHANGE PER ITERATION OF TRBFLX
1665  C    ITRTRB        -         NUMBER OF ITERATIONS OF TRBLFX  C    ITRTRB        -         NUMBER OF ITERATIONS OF TRBFLX
1666  C    KMBG          -         BACKGROUND VALUE OF MOMENTUM TRANSFER COEF  C    KMBG          -         BACKGROUND VALUE OF MOMENTUM TRANSFER COEF
1667  C    KHBG          -         BACKGROUND VALUE OF HEAT TRANSFER COEF  C    KHBG          -         BACKGROUND VALUE OF HEAT TRANSFER COEF
1668  C    NLEV          -         NUMBER OF ATMOSPHERIC LEVELS TO CALCULATE  C    NLEV          -         NUMBER OF ATMOSPHERIC LEVELS TO CALCULATE
# Line 2296  C Line 2476  C
2476        rhokdz(i,nlev) = 0.0              rhokdz(i,nlev) = 0.0      
2477        enddo        enddo
2478    
2479        do nt = 1,ntrace  c     do nt = 1,ntrace
2480        do  i = 1,irun  c     do  i = 1,irun
2481        tracers(i,nlev+1,nt) = tracers(i,nlev,nt)  c     tracers(i,nlev+1,nt) = tracers(i,nlev,nt)
2482        enddo  c     enddo
2483        CALL TRBDIF(tracers(1,1,nt),DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,  c     CALL TRBDIF(tracers(1,1,nt),DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,
2484       .                                         NLEV,4,0. _d 0,irun)  c    .                                         NLEV,4,0. _d 0,irun)
2485        enddo  c     enddo
2486  C  C
2487  C   CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: U AND V  C   CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: U AND V
2488  C  C
# Line 3175  C Local Variables Line 3355  C Local Variables
3355        _RL X0NUM(irun),X1NUM(irun),X0DEN(irun)        _RL X0NUM(irun),X1NUM(irun),X0DEN(irun)
3356        _RL X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)        _RL X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)
3357        _RL cm3,cm4,cm5,cm8        _RL cm3,cm4,cm5,cm8
3358        integer ibit,index        integer ibit,indx
3359        integer i        integer i
3360  C  C
3361        CM3 =   sqrt( 0.2/CM1-0.01 )        CM3 =   sqrt( 0.2/CM1-0.01 )
# Line 3210  C ************************************** Line 3390  C **************************************
3390  C  C
3391        IF(IBIT.LE.0)  GO TO 100        IF(IBIT.LE.0)  GO TO 100
3392  C  C
3393        INDEX = 0        indx = 0
3394        DO 9002 I = 1,IRUN        DO 9002 I = 1,IRUN
3395         IF (INTSTB(I).EQ.1)THEN         IF (INTSTB(I).EQ.1)THEN
3396          INDEX = INDEX + 1          indx = indx + 1
3397          Z(INDEX) = VZZ(I)          Z(indx) = VZZ(I)
3398          Z0(INDEX) = ZZ0(I)          Z0(indx) = ZZ0(I)
3399         ENDIF         ENDIF
3400   9002 CONTINUE   9002 CONTINUE
3401  C  C
# Line 3253  C Line 3433  C
3433         PSI2(I) = PSI2(I) + DX(I)         PSI2(I) = PSI2(I) + DX(I)
3434   9006 CONTINUE   9006 CONTINUE
3435  C  C
3436        INDEX = 0        indx = 0
3437        DO 9008 I = 1,IRUN        DO 9008 I = 1,IRUN
3438         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3439          INDEX = INDEX + 1          indx = indx + 1
3440          VPSIM(I) = PSI2(INDEX)          VPSIM(I) = PSI2(indx)
3441          VX(I) = X1(INDEX)          VX(I) = X1(indx)
3442          VXS(I) = X0(INDEX)          VXS(I) = X0(indx)
3443         ENDIF         ENDIF
3444   9008 CONTINUE   9008 CONTINUE
3445  C  C
# Line 3286  C Line 3466  C
3466         PSI2(I) = PSI2(I) - Y1(I) + Y0(I)         PSI2(I) = PSI2(I) - Y1(I) + Y0(I)
3467   9010 CONTINUE   9010 CONTINUE
3468  C  C
3469        INDEX = 0        indx = 0
3470        DO 9012 I = 1,IRUN        DO 9012 I = 1,IRUN
3471         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3472         INDEX = INDEX + 1         indx = indx + 1
3473         VPSIH(I) = PSI2(INDEX)         VPSIH(I) = PSI2(indx)
3474         VY(I) = Y1(INDEX)         VY(I) = Y1(indx)
3475         VYS(I) = Y0(INDEX)         VYS(I) = Y0(indx)
3476         ENDIF         ENDIF
3477   9012 CONTINUE   9012 CONTINUE
3478  C  C
# Line 3315  C Line 3495  C
3495         ENDIF         ENDIF
3496   9014 CONTINUE   9014 CONTINUE
3497        IF(IBIT.LE.0)  GO TO 300        IF(IBIT.LE.0)  GO TO 300
3498        INDEX = 0        indx = 0
3499  #ifdef CRAY  #ifdef CRAY
3500  CDIR$ NOVECTOR  CDIR$ NOVECTOR
3501  #endif  #endif
3502        DO 9016 I = 1,IRUN        DO 9016 I = 1,IRUN
3503         IF (INTSTB(I).EQ.1)THEN         IF (INTSTB(I).EQ.1)THEN
3504          INDEX = INDEX + 1          indx = indx + 1
3505          Z(INDEX) = VZZ(I)          Z(indx) = VZZ(I)
3506          Z0(INDEX) = ZZ0(I)          Z0(indx) = ZZ0(I)
3507          ARG1(INDEX) = VZH(I)          ARG1(indx) = VZH(I)
3508         ENDIF         ENDIF
3509   9016 CONTINUE   9016 CONTINUE
3510  #ifdef CRAY  #ifdef CRAY
# Line 3381  C Line 3561  C
3561         PSI2(I) = TEMP(I) + CM6 * ARG1(I)         PSI2(I) = TEMP(I) + CM6 * ARG1(I)
3562   9020 CONTINUE   9020 CONTINUE
3563  C  C
3564        INDEX = 0        indx = 0
3565        DO 9030 I = 1,IRUN        DO 9030 I = 1,IRUN
3566         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3567         INDEX = INDEX + 1         indx = indx + 1
3568         VPSIM(I) = PSI2(INDEX)         VPSIM(I) = PSI2(indx)
3569         VX(I) = X1(INDEX)         VX(I) = X1(indx)
3570         VXS(I) = X0(INDEX)         VXS(I) = X0(indx)
3571         ENDIF         ENDIF
3572   9030 CONTINUE   9030 CONTINUE
3573  C  C
# Line 3413  C Line 3593  C
3593         PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)         PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)
3594   9024 CONTINUE   9024 CONTINUE
3595  C  C
3596        INDEX = 0        indx = 0
3597        DO 9026 I = 1,IRUN        DO 9026 I = 1,IRUN
3598         IF( INTSTB(I).EQ.1 ) THEN         IF( INTSTB(I).EQ.1 ) THEN
3599         INDEX = INDEX + 1         indx = indx + 1
3600         VPSIH(I) = PSI2(INDEX)         VPSIH(I) = PSI2(indx)
3601         VY(I) = Y1(INDEX)         VY(I) = Y1(indx)
3602         VYS(I) = X0(INDEX)         VYS(I) = X0(indx)
3603         ENDIF         ENDIF
3604   9026 CONTINUE   9026 CONTINUE
3605  C  C
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22