/[MITgcm]/MITgcm/pkg/kpp/kpp_routines.F
ViewVC logotype

Diff of /MITgcm/pkg/kpp/kpp_routines.F

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

revision 1.7 by cnh, Sun Feb 4 14:38:50 2001 UTC revision 1.22 by heimbach, Thu Feb 10 05:41:44 2005 UTC
# Line 23  c*************************************** Line 23  c***************************************
23       I       mytime, mythid       I       mytime, mythid
24       I     , kmtj, shsq, dvsq, ustar       I     , kmtj, shsq, dvsq, ustar
25       I     , bo, bosol, dbloc, Ritop, coriol       I     , bo, bosol, dbloc, Ritop, coriol
26       I     , ikey       I     , ikppkey
27       O     , diffus       O     , diffus
28       U     , ghat       U     , ghat
29       O     , hbl )       O     , hbl )
# Line 82  c           where hbl(i,j) -> hbl((j-1)* Line 82  c           where hbl(i,j) -> hbl((j-1)*
82        _KPP_RL Ritop (imt,Nr)        _KPP_RL Ritop (imt,Nr)
83        _KPP_RL coriol(imt   )        _KPP_RL coriol(imt   )
84    
85        integer ikey        integer ikppkey
86    
87  c output  c output
88  c     diffus (imt,1)  - vertical viscosity coefficient                  (m^2/s)  c     diffus (imt,1)  - vertical viscosity coefficient                  (m^2/s)
# Line 125  c instability. Line 125  c instability.
125  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
126  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
127    
128  CADJ STORE ghat = comlev1_kpp, key = ikey  cph(
129    cph these storings avoid recomp. of Ri_iwmix
130    CADJ STORE ghat  = comlev1_kpp, key = ikppkey
131    CADJ STORE dbloc = comlev1_kpp, key = ikppkey
132    cph)
133        call Ri_iwmix (        call Ri_iwmix (
134       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
135       I     , ikey       I     , ikppkey
136       O     , diffus )       O     , diffus )
137    
138    cph(
139    cph these storings avoid recomp. of Ri_iwmix
140    cph DESPITE TAFs 'not necessary' warning!
141    CADJ STORE dbloc  = comlev1_kpp, key = ikppkey
142    CADJ STORE shsq   = comlev1_kpp, key = ikppkey
143    CADJ STORE ghat   = comlev1_kpp, key = ikppkey
144    CADJ STORE diffus = comlev1_kpp, key = ikppkey
145    cph)
146    
147  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
148  c set seafloor values to zero and fill extra "Nrp1" coefficients  c set seafloor values to zero and fill extra "Nrp1" coefficients
149  c for blmix  c for blmix
# Line 155  c--------------------------------------- Line 167  c---------------------------------------
167       I       mytime, mythid       I       mytime, mythid
168       I     , kmtj       I     , kmtj
169       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
170       I     , ikey       I     , ikppkey
171       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
172       &     )       &     )
173    
174  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey
175    
176  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
177  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 167  c--------------------------------------- Line 179  c---------------------------------------
179    
180        call blmix (        call blmix (
181       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
182       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
183       &     )       &     )
184    cph(
185  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey
186    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey
187    cph)
188    
189  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
190  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 181  c--------------------------------------- Line 195  c---------------------------------------
195       U     , ghat       U     , ghat
196       O     , blmc )       O     , blmc )
197    
198    cph(
199    cph avoids recomp. of enhance
200    CADJ STORE blmc = comlev1_kpp, key = ikppkey
201    cph)
202    
203  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
204  c combine interior and boundary layer coefficients and nonlocal term  c combine interior and boundary layer coefficients and nonlocal term
205    c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
206    c (< 1 level), diffusivity blmc can become negative.  The max-s below
207    c are a hack until this problem is properly diagnosed and fixed.
208  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
   
209        do k = 1, Nr        do k = 1, Nr
210           do i = 1, imt           do i = 1, imt
211              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
212                 do md = 1, mdiff                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
213                    diffus(i,k,md) = blmc(i,k,md)                 diffus(i,k,2) = max ( blmc(i,k,2), diffKrNrS(k) )
214                 end do                 diffus(i,k,3) = max ( blmc(i,k,3), diffKrNrT(k) )
215              else              else
216                 ghat(i,k) = 0.                 ghat(i,k) = 0.
217              endif              endif
# Line 208  c*************************************** Line 229  c***************************************
229       I       mytime, mythid       I       mytime, mythid
230       I     , kmtj       I     , kmtj
231       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
232       I     , ikey       I     , ikppkey
233       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
234       &     )       &     )
235    
# Line 265  c coriol    : Coriolis parameter Line 286  c coriol    : Coriolis parameter
286        _KPP_RL bo    (imt)        _KPP_RL bo    (imt)
287        _KPP_RL bosol (imt)        _KPP_RL bosol (imt)
288        _KPP_RL coriol(imt)        _KPP_RL coriol(imt)
289        integer ikey        integer ikppkey, kkppkey
290    
291  c  output  c  output
292  c--------  c--------
# Line 319  c     initialize hbl and kbl to bottomed Line 340  c     initialize hbl and kbl to bottomed
340    
341        do kl = 2, Nr        do kl = 2, Nr
342    
343    #ifdef ALLOW_AUTODIFF_TAMC
344             kkppkey = (ikppkey-1)*Nr + kl
345    #endif
346    
347  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
348    
349           do i = 1, imt           do i = 1, imt
350              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
351           end do           end do
352    CADJ store worka = comlev1_kpp_k, key = kkppkey
353           call SWFRAC(           call SWFRAC(
354       I        imt, hbf,       I        imt, hbf,
355       I        mytime, mythid,       I        mytime, mythid,
356       U        worka )       U        worka )
357    CADJ store worka = comlev1_kpp_k, key = kkppkey
358    
359    
360           do i = 1, imt           do i = 1, imt
361    
# Line 350  c--------------------------------------- Line 378  c---------------------------------------
378           call wscale (           call wscale (
379       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
380       O        wm, ws )       O        wm, ws )
381    CADJ store ws = comlev1_kpp_k, key = kkppkey
382    
383           do i = 1, imt           do i = 1, imt
384    
# Line 385  c Line 414  c
414    
415           end do           end do
416        end do        end do
417          
418    cph(
419    cph  without this store, there is a recomputation error for
420    cph  rib in adbldepth (probably partial recomputation problem)    
421    CADJ store Rib = comlev1_kpp
422    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
423    cph)
424    
425        do kl = 2, Nr        do kl = 2, Nr
426           do i = 1, imt           do i = 1, imt
427              if (kbl(i).eq.kmtj(i) .and. Rib(i,kl).gt.Ricr) kbl(i) = kl              if (kbl(i).eq.kmtj(i) .and. Rib(i,kl).gt.Ricr) kbl(i) = kl
# Line 393  c Line 429  c
429        end do        end do
430    
431  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
432  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
433    
434        do i = 1, imt        do i = 1, imt
435           kl = kbl(i)           kl = kbl(i)
# Line 406  c     linearly interpolate to find hbl w Line 442  c     linearly interpolate to find hbl w
442        end do        end do
443    
444  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
445  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
446    
447  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
448  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
# Line 415  c--------------------------------------- Line 451  c---------------------------------------
451        do i = 1, imt        do i = 1, imt
452           worka(i) = hbl(i)           worka(i) = hbl(i)
453        end do        end do
454    CADJ store worka = comlev1_kpp
455    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
456        call SWFRAC(        call SWFRAC(
457       I     imt, minusone,       I     imt, minusone,
458       I     mytime, mythid,       I     mytime, mythid,
459       U     worka )       U     worka )
460    CADJ store worka = comlev1_kpp
461    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
462    
463        do i = 1, imt        do i = 1, imt
464           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
465        end do        end do
466  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
467  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
468    
469  c--   ensure bfsfc is never 0  c--   ensure bfsfc is never 0
470        do i = 1, imt        do i = 1, imt
# Line 432  c--   ensure bfsfc is never 0 Line 472  c--   ensure bfsfc is never 0
472           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
473        end do        end do
474    
475  CADJ store bfsfc = comlev1_kpp  cph(
476  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  cph  added stable to store list to avoid extensive recomp.
477    CADJ store bfsfc, stable = comlev1_kpp
478    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
479    cph)
480    
481  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
482  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
# Line 451  c--------------------------------------- Line 494  c---------------------------------------
494           end if           end if
495        end do        end do
496  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
497  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
498    
499        do i = 1, imt        do i = 1, imt
500           hbl(i) = max(hbl(i),minKPPhbl)           hbl(i) = max(hbl(i),minKPPhbl)
# Line 459  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 502  CADJ &   , key = ikey, shape = (/ (sNx+2
502        end do        end do
503    
504  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
505  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
506    
507  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
508  c      find new kbl  c      find new kbl
# Line 480  c--------------------------------------- Line 523  c---------------------------------------
523        do i = 1, imt        do i = 1, imt
524           worka(i) = hbl(i)           worka(i) = hbl(i)
525        end do        end do
526    CADJ store worka = comlev1_kpp
527    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
528        call SWFRAC(        call SWFRAC(
529       I     imt, minusone,       I     imt, minusone,
530       I     mytime, mythid,       I     mytime, mythid,
531       U     worka )       U     worka )
532    CADJ store worka = comlev1_kpp
533    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
534    
535        do i = 1, imt        do i = 1, imt
536           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
537        end do        end do
538  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
539  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
540    
541  c--   ensures bfsfc is never 0  c--   ensures bfsfc is never 0
542        do i = 1, imt        do i = 1, imt
# Line 611  c*************************************** Line 658  c***************************************
658    
659        subroutine Ri_iwmix (        subroutine Ri_iwmix (
660       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm
661       I     , ikey       I     , ikppkey
662       O     , diffus )       O     , diffus )
663    
664  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
# Line 635  c     dblocSm(imt,Nr)      horizontally Line 682  c     dblocSm(imt,Nr)      horizontally
682        _KPP_RL shsq   (imt,Nr)        _KPP_RL shsq   (imt,Nr)
683        _KPP_RL dbloc  (imt,Nr)        _KPP_RL dbloc  (imt,Nr)
684        _KPP_RL dblocSm(imt,Nr)        _KPP_RL dblocSm(imt,Nr)
685        integer ikey        integer ikppkey
686    
687  c output  c output
688  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
# Line 651  c     fRi, fcon             function of Line 698  c     fRi, fcon             function of
698        _KPP_RL Rig        _KPP_RL Rig
699        _KPP_RL fRi, fcon        _KPP_RL fRi, fcon
700        _KPP_RL ratio        _KPP_RL ratio
701        integer i, ki, mr        integer i, ki
702        _KPP_RL c1, c0        _KPP_RL c1, c0
703    
704  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
705          integer mr
706  CADJ INIT kpp_ri_tape_mr = common, 1  CADJ INIT kpp_ri_tape_mr = common, 1
707  #endif  #endif
708    
# Line 670  c     set values at bottom and below to Line 718  c     set values at bottom and below to
718  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
719  C     break data flow dependence on diffus  C     break data flow dependence on diffus
720        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
721    
722          do ki = 1, Nr
723             do i = 1, imt
724                diffus(i,ki,1) = 0.
725                diffus(i,ki,2) = 0.
726                diffus(i,ki,3) = 0.
727             enddo
728          enddo
729  #endif  #endif
730                
731    
732        do ki = 1, Nr        do ki = 1, Nr
733           do i = 1, imt           do i = 1, imt
734              if     (kmtj(i) .EQ. 0      ) then              if     (kmtj(i) .LE. 1      ) then
735                 diffus(i,ki,1) = 0.                 diffus(i,ki,1) = 0.
736                 diffus(i,ki,2) = 0.                 diffus(i,ki,2) = 0.
737              elseif (ki      .GE. kmtj(i)) then              elseif (ki      .GE. kmtj(i)) then
# Line 687  C     break data flow dependence on diff Line 744  C     break data flow dependence on diff
744              endif              endif
745           end do           end do
746        end do        end do
747    CADJ store diffus = comlev1_kpp, key = ikppkey
748    
749  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
750  c     vertically smooth Ri  c     vertically smooth Ri
# Line 724  c  evaluate f of smooth Ri for shear ins Line 782  c  evaluate f of smooth Ri for shear ins
782  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
783  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
784  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
785    
786              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0  #ifndef EXCLUDE_KPP_SHEAR_MIX
787              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0              if ( .NOT. inAdMode ) then
788              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
789                   diffus(i,ki,2) = diffKrNrS(ki)+ fcon*difscon + fRi*difs0
790                   diffus(i,ki,3) = diffKrNrT(ki)+ fcon*difscon + fRi*difs0
791                else
792                   diffus(i,ki,1) = viscAr
793                   diffus(i,ki,2) = diffKrNrS(ki)
794                   diffus(i,ki,3) = diffKrNrT(ki)
795                endif
796    #else
797                diffus(i,ki,1) = viscAr
798                diffus(i,ki,2) = diffKrNrS(ki)
799                diffus(i,ki,3) = diffKrNrT(ki)
800    #endif
801    
802           end do           end do
803        end do        end do
# Line 838  c     Apply horizontal smoothing to KPP Line 908  c     Apply horizontal smoothing to KPP
908    
909        IMPLICIT NONE        IMPLICIT NONE
910  #include "SIZE.h"  #include "SIZE.h"
911    #include "GRID.h"
912  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
913    
914  c     input  c     input
# Line 869  c     local Line 940  c     local
940              im1 = i-1              im1 = i-1
941              ip1 = i+1              ip1 = i+1
942              tempVar =              tempVar =
943       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
944       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
945       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
946       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
947       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
948       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
949       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
950       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
951       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
952              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
953                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
954       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
955       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
956       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
957       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
958       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
959       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
960       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
961       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
962       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
963       &              / tempVar       &              / tempVar
964              ELSE              ELSE
965                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 918  c     Apply horizontal smoothing to glob Line 989  c     Apply horizontal smoothing to glob
989    
990        IMPLICIT NONE        IMPLICIT NONE
991  #include "SIZE.h"  #include "SIZE.h"
992    #include "GRID.h"
993  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
994    
995  c     input  c     input
# Line 949  c     local Line 1021  c     local
1021              im1 = i-1              im1 = i-1
1022              ip1 = i+1              ip1 = i+1
1023              tempVar =              tempVar =
1024       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
1025       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
1026       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
1027       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
1028       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
1029       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
1030       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
1031       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
1032       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
1033              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
1034                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
1035       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
1036       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
1037       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
1038       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
1039       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
1040       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
1041       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
1042       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
1043       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
1044       &              / tempVar       &              / tempVar
1045              ELSE              ELSE
1046                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 992  c*************************************** Line 1064  c***************************************
1064    
1065        subroutine blmix (        subroutine blmix (
1066       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1067       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
1068       &     )       &     )
1069    
1070  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
# Line 1033  c     sigma(imt)                  normal Line 1105  c     sigma(imt)                  normal
1105        _KPP_RL blmc (imt,Nr,mdiff)        _KPP_RL blmc (imt,Nr,mdiff)
1106        _KPP_RL ghat (imt,Nr)        _KPP_RL ghat (imt,Nr)
1107        _KPP_RL sigma(imt)        _KPP_RL sigma(imt)
1108        integer ikey        integer ikppkey, kkppkey
1109    
1110  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1111    
# Line 1062  c--------------------------------------- Line 1134  c---------------------------------------
1134           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1135        end do        end do
1136    
1137    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1138        call wscale (        call wscale (
1139       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1140       O        wm, ws )       O        wm, ws )
1141    CADJ STORE wm = comlev1_kpp, key = ikppkey
1142    CADJ STORE ws = comlev1_kpp, key = ikppkey
1143    
1144        do i = 1, imt        do i = 1, imt
1145           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1146           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1147        end do        end do
1148  CADJ STORE wm = comlev1_kpp, key = ikey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1149  CADJ STORE ws = comlev1_kpp, key = ikey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1150    
1151        do i = 1, imt        do i = 1, imt
1152    
# Line 1107  c--------------------------------------- Line 1182  c---------------------------------------
1182       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1183           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1184           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
          dat1m(i) = min(dat1m(i),p0)  
1185    
1186           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1187           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
          dat1s(i) = min(dat1s(i),p0)  
1188    
1189           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1190           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
          dat1t(i) = min(dat1t(i),p0)  
1191    
1192        end do        end do
1193    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1194    CADJ STORE gat1m = comlev1_kpp, key = ikppkey
1195    CADJ STORE gat1s = comlev1_kpp, key = ikppkey
1196    CADJ STORE gat1t = comlev1_kpp, key = ikppkey
1197    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1198    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1199    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1200    #endif
1201          do i = 1, imt
1202             dat1m(i) = min(dat1m(i),p0)
1203             dat1s(i) = min(dat1s(i),p0)
1204             dat1t(i) = min(dat1t(i),p0)
1205          end do
1206    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1207    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1208    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1209    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1210    #endif
1211    
1212        do ki = 1, Nr        do ki = 1, Nr
1213    
1214    #ifdef ALLOW_AUTODIFF_TAMC
1215             kkppkey = (ikppkey-1)*Nr + ki
1216    #endif
1217    
1218  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1219  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1220  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1129  c--------------------------------------- Line 1223  c---------------------------------------
1223              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1224              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1225           end do           end do
1226    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1227    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1228    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1229    #endif
1230    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1231           call wscale (           call wscale (
1232       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1233       O        wm, ws )       O        wm, ws )
1234    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1235    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1236    
1237  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1238  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1175  c--------------------------------------- Line 1276  c---------------------------------------
1276       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1277        end do        end do
1278    
1279    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1280    CADJ STORE wm = comlev1_kpp, key = ikppkey
1281    CADJ STORE ws = comlev1_kpp, key = ikppkey
1282    #endif
1283    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1284        call wscale (        call wscale (
1285       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1286       O        wm, ws )       O        wm, ws )
1287    CADJ STORE wm = comlev1_kpp, key = ikppkey
1288    CADJ STORE ws = comlev1_kpp, key = ikppkey
1289    
1290        do i = 1, imt        do i = 1, imt
1291           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1264  c     fraction hbl lies beteen zgrid nei Line 1372  c     fraction hbl lies beteen zgrid nei
1372  c*************************************************************************  c*************************************************************************
1373    
1374        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1375       I     bi, bj, myThid,       I     ikppkey, bi, bj, myThid,
1376       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)
1377  c  c
1378  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1326  c     work1, work2 - work arrays for hol Line 1434  c     work1, work2 - work arrays for hol
1434        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1435        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1436        INTEGER I, J, K        INTEGER I, J, K
1437          INTEGER ikppkey, kkppkey
1438    
1439  c calculate density, alpha, beta in surface layer, and set dbsfc to zero  c calculate density, alpha, beta in surface layer, and set dbsfc to zero
1440    
1441          kkppkey = (ikppkey-1)*Nr + 1
1442    
1443    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1444    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1445    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1446    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1447        call FIND_RHO(        call FIND_RHO(
1448       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,
1449       I     theta, salt,       I     theta, salt,
1450       O     WORK1,       O     WORK1,
1451       I     myThid )       I     myThid )
1452    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1453    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1454    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1455    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1456    
1457        call FIND_ALPHA(        call FIND_ALPHA(
1458       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,
1459       O     WORK2 )       O     WORK2 )
1460    
1461        call FIND_BETA(        call FIND_BETA(
1462       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,
1463       O     WORK3 )       O     WORK3 )
1464    
1465        DO J = jbot, jtop        DO J = jbot, jtop
1466           DO I = ibot, itop           DO I = ibot, itop
1467              RHO1(I,J)      = WORK1(I,J) + rhonil              RHO1(I,J)      = WORK1(I,J) + rhoConst
1468              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1469              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
1470              DBSFC(I,J,1)   = 0.              DBSFC(I,J,1)   = 0.
# Line 1357  c calculate alpha, beta, and gradients i Line 1476  c calculate alpha, beta, and gradients i
1476  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1477        DO K = 2, Nr        DO K = 2, Nr
1478    
1479          kkppkey = (ikppkey-1)*Nr + k
1480    
1481    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1482    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1483    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1484    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1485           call FIND_RHO(           call FIND_RHO(
1486       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, ibot, itop, jbot, jtop, K, K,
1487       I        theta, salt,       I        theta, salt,
1488       O        RHOK,       O        RHOK,
1489       I        myThid )       I        myThid )
1490    
1491    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1492    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1493    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1494    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1495           call FIND_RHO(           call FIND_RHO(
1496       I        bi, bj, ibot, itop, jbot, jtop, K-1, K, eosType,       I        bi, bj, ibot, itop, jbot, jtop, K-1, K,
1497       I        theta, salt,       I        theta, salt,
1498       O        RHOKM1,       O        RHOKM1,
1499       I        myThid )       I        myThid )
1500    
1501    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1502    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1503    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1504    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1505           call FIND_RHO(           call FIND_RHO(
1506       I        bi, bj, ibot, itop, jbot, jtop, 1, K, eosType,       I        bi, bj, ibot, itop, jbot, jtop, 1, K,
1507       I        theta, salt,       I        theta, salt,
1508       O        RHO1K,       O        RHO1K,
1509       I        myThid )       I        myThid )
1510    
1511    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1512    CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey
1513    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k, key=kkppkey
1514    CADJ STORE rho1k (:,:)          = comlev1_kpp_k, key=kkppkey
1515    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1516    
1517           call FIND_ALPHA(           call FIND_ALPHA(
1518       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, ibot, itop, jbot, jtop, K, K,
1519       O        WORK1 )       O        WORK1 )
1520    
1521           call FIND_BETA(           call FIND_BETA(
1522       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, ibot, itop, jbot, jtop, K, K,
1523       O        WORK2 )       O        WORK2 )
1524    
1525           DO J = jbot, jtop           DO J = jbot, jtop
# Line 1388  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO Line 1527  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1527                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1528                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1529                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
1530       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1531                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /
1532       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1533              END DO              END DO
1534           END DO           END DO
1535    
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22