/[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.18 by heimbach, Tue Dec 16 20:58:57 2003 UTC revision 1.26 by dimitri, Mon Apr 23 20:46:49 2007 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     , diffusKzS, diffusKzT
27       I     , ikppkey       I     , ikppkey
28       O     , diffus       O     , diffus
29       U     , ghat       U     , ghat
# Line 47  c--------------------------------------- Line 48  c---------------------------------------
48  #include "SIZE.h"  #include "SIZE.h"
49  #include "EEPARAMS.h"  #include "EEPARAMS.h"
50  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
51  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
52    
53  c input  c input
54  c     myTime          - current time in simulation  c     myTime           - current time in simulation
55  c     myThid          - thread number for this instance of the routine  c     myThid           - thread number for this instance of the routine
56  c     kmtj   (imt)    - number of vertical layers on this row  c     kmtj   (imt)     - number of vertical layers on this row
57  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
58  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
59  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     ustar  (imt)     - surface friction velocity                        (m/s)
60  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
61  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
62  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
63  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
64  c                         stored in ghat to save space  c                          stored in ghat to save space
65  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
66  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
67  c     coriol (imt)    - Coriolis parameter                                (1/s)  c     coriol (imt)     - Coriolis parameter                               (1/s)
68    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
69    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
70  c     note: there is a conversion from 2-D to 1-D for input output variables,  c     note: there is a conversion from 2-D to 1-D for input output variables,
71  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
72  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
73    
74        _RL     mytime        _RL     mytime
75        integer mythid        integer mythid
76        integer kmtj  (imt   )        integer kmtj     (imt   )
77        _KPP_RL shsq  (imt,Nr)        _KPP_RL shsq     (imt,Nr)
78        _KPP_RL dvsq  (imt,Nr)        _KPP_RL dvsq     (imt,Nr)
79        _KPP_RL ustar (imt   )        _KPP_RL ustar    (imt   )
80        _KPP_RL bo    (imt   )        _KPP_RL bo       (imt   )
81        _KPP_RL bosol (imt   )        _KPP_RL bosol    (imt   )
82        _KPP_RL dbloc (imt,Nr)        _KPP_RL dbloc    (imt,Nr)
83        _KPP_RL Ritop (imt,Nr)        _KPP_RL Ritop    (imt,Nr)
84        _KPP_RL coriol(imt   )        _KPP_RL coriol   (imt   )
85          _RL     diffusKzS(imt,Nr)
86          _RL     diffusKzT(imt,Nr)
87    
88        integer ikppkey        integer ikppkey
89    
# Line 132  CADJ STORE dbloc = comlev1_kpp, key = ik Line 135  CADJ STORE dbloc = comlev1_kpp, key = ik
135  cph)  cph)
136        call Ri_iwmix (        call Ri_iwmix (
137       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
138         I     , diffusKzS, diffusKzT
139       I     , ikppkey       I     , ikppkey
140       O     , diffus )       O     , diffus )
141    
# Line 210  c--------------------------------------- Line 214  c---------------------------------------
214           do i = 1, imt           do i = 1, imt
215              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
216                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
217                 diffus(i,k,2) = max ( blmc(i,k,2), diffKrS )                 diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
218                 diffus(i,k,3) = max ( blmc(i,k,3), diffKrT )                 diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
219              else              else
220                 ghat(i,k) = 0.                 ghat(i,k) = 0.
221              endif              endif
# Line 257  c Line 261  c
261        IMPLICIT NONE        IMPLICIT NONE
262    
263  #include "SIZE.h"  #include "SIZE.h"
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
264  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
 #include "FFIELDS.h"  
265    
266  c input  c input
267  c------  c------
# Line 286  c coriol    : Coriolis parameter Line 287  c coriol    : Coriolis parameter
287        _KPP_RL bo    (imt)        _KPP_RL bo    (imt)
288        _KPP_RL bosol (imt)        _KPP_RL bosol (imt)
289        _KPP_RL coriol(imt)        _KPP_RL coriol(imt)
290        integer ikppkey        integer ikppkey, kkppkey
291    
292  c  output  c  output
293  c--------  c--------
# Line 340  c     initialize hbl and kbl to bottomed Line 341  c     initialize hbl and kbl to bottomed
341    
342        do kl = 2, Nr        do kl = 2, Nr
343    
344    #ifdef ALLOW_AUTODIFF_TAMC
345             kkppkey = (ikppkey-1)*Nr + kl
346    #endif
347    
348  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
349    
350           do i = 1, imt           do i = 1, imt
351              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
352           end do           end do
353    CADJ store worka = comlev1_kpp_k, key = kkppkey
354           call SWFRAC(           call SWFRAC(
355       I        imt, hbf,       I        imt, hbf,
356       I        mytime, mythid,       I        mytime, mythid,
357       U        worka )       U        worka )
358    CADJ store worka = comlev1_kpp_k, key = kkppkey
359    
360    
361           do i = 1, imt           do i = 1, imt
362    
# Line 371  c--------------------------------------- Line 379  c---------------------------------------
379           call wscale (           call wscale (
380       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
381       O        wm, ws )       O        wm, ws )
382    CADJ store ws = comlev1_kpp_k, key = kkppkey
383    
384           do i = 1, imt           do i = 1, imt
385    
# Line 443  c--------------------------------------- Line 452  c---------------------------------------
452        do i = 1, imt        do i = 1, imt
453           worka(i) = hbl(i)           worka(i) = hbl(i)
454        end do        end do
455    CADJ store worka = comlev1_kpp
456    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
457        call SWFRAC(        call SWFRAC(
458       I     imt, minusone,       I     imt, minusone,
459       I     mytime, mythid,       I     mytime, mythid,
460       U     worka )       U     worka )
461    CADJ store worka = comlev1_kpp
462    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
463    
464        do i = 1, imt        do i = 1, imt
465           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
# Line 511  c--------------------------------------- Line 524  c---------------------------------------
524        do i = 1, imt        do i = 1, imt
525           worka(i) = hbl(i)           worka(i) = hbl(i)
526        end do        end do
527    CADJ store worka = comlev1_kpp
528    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
529        call SWFRAC(        call SWFRAC(
530       I     imt, minusone,       I     imt, minusone,
531       I     mytime, mythid,       I     mytime, mythid,
532       U     worka )       U     worka )
533    CADJ store worka = comlev1_kpp
534    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
535    
536        do i = 1, imt        do i = 1, imt
537           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
# Line 642  c*************************************** Line 659  c***************************************
659    
660        subroutine Ri_iwmix (        subroutine Ri_iwmix (
661       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm
662         I     , diffusKzS, diffusKzT
663       I     , ikppkey       I     , ikppkey
664       O     , diffus )       O     , diffus )
665    
# Line 662  c     kmtj   (imt)         number of ver Line 680  c     kmtj   (imt)         number of ver
680  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
681  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
682  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
683    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
684    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
685        integer kmtj   (imt)        integer kmtj   (imt)
686        _KPP_RL shsq   (imt,Nr)        _KPP_RL shsq   (imt,Nr)
687        _KPP_RL dbloc  (imt,Nr)        _KPP_RL dbloc  (imt,Nr)
688        _KPP_RL dblocSm(imt,Nr)        _KPP_RL dblocSm(imt,Nr)
689          _RL diffusKzS(imt,Nr)
690          _RL diffusKzT(imt,Nr)
691        integer ikppkey        integer ikppkey
692    
693  c output  c output
# Line 702  c     set values at bottom and below to Line 724  c     set values at bottom and below to
724  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
725  C     break data flow dependence on diffus  C     break data flow dependence on diffus
726        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
727    
728          do ki = 1, Nr
729             do i = 1, imt
730                diffus(i,ki,1) = 0.
731                diffus(i,ki,2) = 0.
732                diffus(i,ki,3) = 0.
733             enddo
734          enddo
735  #endif  #endif
736                
737    
738        do ki = 1, Nr        do ki = 1, Nr
739           do i = 1, imt           do i = 1, imt
# Line 719  C     break data flow dependence on diff Line 750  C     break data flow dependence on diff
750              endif              endif
751           end do           end do
752        end do        end do
753    CADJ store diffus = comlev1_kpp, key = ikppkey
754    
755  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
756  c     vertically smooth Ri  c     vertically smooth Ri
# Line 758  c            evaluate diffusivities and Line 790  c            evaluate diffusivities and
790  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
791    
792  #ifndef EXCLUDE_KPP_SHEAR_MIX  #ifndef EXCLUDE_KPP_SHEAR_MIX
793              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0              if ( .NOT. inAdMode ) then
794              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
795              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0                 diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
796                   diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
797                else
798                   diffus(i,ki,1) = viscAr
799                   diffus(i,ki,2) = diffusKzS(i,ki)
800                   diffus(i,ki,3) = diffusKzT(i,ki)
801                endif
802  #else  #else
803              diffus(i,ki,1) = viscAr              diffus(i,ki,1) = viscAr
804              diffus(i,ki,2) = diffKrS              diffus(i,ki,2) = diffusKzS(i,ki)
805              diffus(i,ki,3) = diffKrT              diffus(i,ki,3) = diffusKzT(i,ki)
806  #endif  #endif
807    
808           end do           end do
# Line 876  c     Apply horizontal smoothing to KPP Line 914  c     Apply horizontal smoothing to KPP
914    
915        IMPLICIT NONE        IMPLICIT NONE
916  #include "SIZE.h"  #include "SIZE.h"
917    #include "GRID.h"
918  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
919    
920  c     input  c     input
# Line 885  c     k      : vertical index used for m Line 924  c     k      : vertical index used for m
924    
925  c     input/output  c     input/output
926  c     fld    : 2-D array to be smoothed  c     fld    : 2-D array to be smoothed
927        _KPP_RL fld( ibot:itop, jbot:jtop )        _KPP_RL fld( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
928    
929  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
930    
931  c     local  c     local
932        integer i, j, im1, ip1, jm1, jp1        integer i, j, im1, ip1, jm1, jp1
933        _KPP_RL tempVar        _KPP_RL tempVar
934        _KPP_RL fld_tmp( ibot:itop, jbot:jtop )        _KPP_RL fld_tmp( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
935    
936        integer    imin       , imax       , jmin       , jmax        integer   imin      , imax          , jmin      , jmax
937        parameter( imin=ibot+1, imax=itop-1, jmin=jbot+1, jmax=jtop-1 )        parameter(imin=2-OLx, imax=sNx+OLx-1, jmin=2-OLy, jmax=sNy+OLy-1)
938    
939        _KPP_RL    p0    , p5    , p25     , p125      , p0625        _KPP_RL    p0    , p5    , p25     , p125      , p0625
940        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
# Line 907  c     local Line 946  c     local
946              im1 = i-1              im1 = i-1
947              ip1 = i+1              ip1 = i+1
948              tempVar =              tempVar =
949       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
950       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
951       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
952       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
953       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
954       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
955       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
956       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
957       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
958              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
959                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
960       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
961       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
962       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
963       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
964       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
965       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
966       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
967       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
968       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
969       &              / tempVar       &              / tempVar
970              ELSE              ELSE
971                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 956  c     Apply horizontal smoothing to glob Line 995  c     Apply horizontal smoothing to glob
995    
996        IMPLICIT NONE        IMPLICIT NONE
997  #include "SIZE.h"  #include "SIZE.h"
998    #include "GRID.h"
999  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1000    
1001  c     input  c     input
# Line 987  c     local Line 1027  c     local
1027              im1 = i-1              im1 = i-1
1028              ip1 = i+1              ip1 = i+1
1029              tempVar =              tempVar =
1030       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
1031       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
1032       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
1033       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
1034       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
1035       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
1036       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
1037       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
1038       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
1039              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
1040                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
1041       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
1042       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
1043       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
1044       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
1045       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
1046       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
1047       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
1048       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
1049       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
1050       &              / tempVar       &              / tempVar
1051              ELSE              ELSE
1052                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 1071  c     sigma(imt)                  normal Line 1111  c     sigma(imt)                  normal
1111        _KPP_RL blmc (imt,Nr,mdiff)        _KPP_RL blmc (imt,Nr,mdiff)
1112        _KPP_RL ghat (imt,Nr)        _KPP_RL ghat (imt,Nr)
1113        _KPP_RL sigma(imt)        _KPP_RL sigma(imt)
1114        integer ikppkey        integer ikppkey, kkppkey
1115    
1116  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1117    
# Line 1100  c--------------------------------------- Line 1140  c---------------------------------------
1140           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1141        end do        end do
1142    
1143    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1144        call wscale (        call wscale (
1145       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1146       O        wm, ws )       O        wm, ws )
1147    CADJ STORE wm = comlev1_kpp, key = ikppkey
1148    CADJ STORE ws = comlev1_kpp, key = ikppkey
1149    
1150        do i = 1, imt        do i = 1, imt
1151           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
# Line 1145  c--------------------------------------- Line 1188  c---------------------------------------
1188       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1189           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1190           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
          dat1m(i) = min(dat1m(i),p0)  
1191    
1192           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1193           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
          dat1s(i) = min(dat1s(i),p0)  
1194    
1195           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1196           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
          dat1t(i) = min(dat1t(i),p0)  
1197    
1198        end do        end do
1199    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1200    CADJ STORE gat1m = comlev1_kpp, key = ikppkey
1201    CADJ STORE gat1s = comlev1_kpp, key = ikppkey
1202    CADJ STORE gat1t = comlev1_kpp, key = ikppkey
1203    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1204    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1205    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1206    #endif
1207          do i = 1, imt
1208             dat1m(i) = min(dat1m(i),p0)
1209             dat1s(i) = min(dat1s(i),p0)
1210             dat1t(i) = min(dat1t(i),p0)
1211          end do
1212    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1213    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1214    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1215    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1216    #endif
1217    
1218        do ki = 1, Nr        do ki = 1, Nr
1219    
1220    #ifdef ALLOW_AUTODIFF_TAMC
1221             kkppkey = (ikppkey-1)*Nr + ki
1222    #endif
1223    
1224  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1225  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1226  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1167  c--------------------------------------- Line 1229  c---------------------------------------
1229              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1230              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1231           end do           end do
1232    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1233    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1234    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1235    #endif
1236    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1237           call wscale (           call wscale (
1238       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1239       O        wm, ws )       O        wm, ws )
1240    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1241    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1242    
1243  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1244  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1213  c--------------------------------------- Line 1282  c---------------------------------------
1282       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1283        end do        end do
1284    
1285    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1286    CADJ STORE wm = comlev1_kpp, key = ikppkey
1287    CADJ STORE ws = comlev1_kpp, key = ikppkey
1288    #endif
1289    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1290        call wscale (        call wscale (
1291       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1292       O        wm, ws )       O        wm, ws )
1293    CADJ STORE wm = comlev1_kpp, key = ikppkey
1294    CADJ STORE ws = comlev1_kpp, key = ikppkey
1295    
1296        do i = 1, imt        do i = 1, imt
1297           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1302  c     fraction hbl lies beteen zgrid nei Line 1378  c     fraction hbl lies beteen zgrid nei
1378  c*************************************************************************  c*************************************************************************
1379    
1380        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1381       I     bi, bj, myThid,       I     ikppkey, bi, bj, myThid,
1382       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)
1383  c  c
1384  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1328  c Line 1404  c
1404  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1405  c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV  c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV
1406  c  c
1407    c     28 april 05: added computation of mixed layer depth KPPmld
1408    c     for a deltaRho equivalent to a deltaTheta of 0.8 deg C
1409    
1410  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1411    
1412        IMPLICIT NONE        IMPLICIT NONE
# Line 1337  c--------------------------------------- Line 1416  c---------------------------------------
1416  #include "PARAMS.h"  #include "PARAMS.h"
1417  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1418  #include "DYNVARS.h"  #include "DYNVARS.h"
1419    #include "GRID.h"
1420    
1421  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1422        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1423        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _KPP_RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1424        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _KPP_RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1425        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _KPP_RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1426        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _KPP_RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1427        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _KPP_RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1428    
1429  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1430    
# Line 1355  c Line 1435  c
1435  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1436  c     rhokm1       - density of t(k-1) & s(k-1) at depth k  c     rhokm1       - density of t(k-1) & s(k-1) at depth k
1437  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1438  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1439    
1440        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1441        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 1363  c     work1, work2 - work arrays for hol Line 1443  c     work1, work2 - work arrays for hol
1443        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1444        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1445        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1446    #ifdef ALLOW_DIAGNOSTICS
1447    c     KPPMLD - mixed layer depth based on density criterion
1448          _RL KPPMLD(1    :sNx    ,1    :sNy    )
1449    #endif /* ALLOW_DIAGNOSTICS */
1450    
1451        INTEGER I, J, K        INTEGER I, J, K
1452          INTEGER ikppkey, kkppkey
1453    
1454  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
1455    
1456          kkppkey = (ikppkey-1)*Nr + 1
1457    
1458    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1459    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1460    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1461    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1462        call FIND_RHO(        call FIND_RHO(
1463       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1464       I     theta, salt,       I     theta, salt,
1465       O     WORK1,       O     WORK1,
1466       I     myThid )       I     myThid )
1467    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1468    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1469    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1470    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1471    
1472        call FIND_ALPHA(        call FIND_ALPHA(
1473       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1474       O     WORK2 )       O     WORK2 )
1475    
1476        call FIND_BETA(        call FIND_BETA(
1477       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1478       O     WORK3 )       O     WORK3 )
1479    
1480        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1481           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1482              RHO1(I,J)      = WORK1(I,J) + rhoConst              RHO1(I,J)      = WORK1(I,J) + rhoConst
1483              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1484              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
# Line 1390  c calculate density, alpha, beta in surf Line 1486  c calculate density, alpha, beta in surf
1486           END DO           END DO
1487        END DO        END DO
1488    
1489    #ifdef ALLOW_DIAGNOSTICS
1490    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1491          IF ( useDiagnostics ) THEN
1492             DO J = 1, sNy
1493                DO I = 1, sNx
1494                   KPPMLD(I,J) = ABS(R_low(I,J,bi,bj))
1495                   WORK3 (I,J) = WORK1(I,J) - 0.8 * WORK2(I,J)
1496                END DO
1497             END DO
1498          ENDIF
1499    #endif /* ALLOW_DIAGNOSTICS */
1500    
1501  c calculate alpha, beta, and gradients in interior layers  c calculate alpha, beta, and gradients in interior layers
1502    
1503  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1504        DO K = 2, Nr        DO K = 2, Nr
1505    
1506          kkppkey = (ikppkey-1)*Nr + k
1507    
1508    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1509    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1510    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1511    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1512           call FIND_RHO(           call FIND_RHO(
1513       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1514       I        theta, salt,       I        theta, salt,
1515       O        RHOK,       O        RHOK,
1516       I        myThid )       I        myThid )
1517    
1518    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1519    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1520    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1521    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1522           call FIND_RHO(           call FIND_RHO(
1523       I        bi, bj, ibot, itop, jbot, jtop, K-1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K-1, K,
1524       I        theta, salt,       I        theta, salt,
1525       O        RHOKM1,       O        RHOKM1,
1526       I        myThid )       I        myThid )
1527    
1528    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1529    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1530    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1531    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1532           call FIND_RHO(           call FIND_RHO(
1533       I        bi, bj, ibot, itop, jbot, jtop, 1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K,
1534       I        theta, salt,       I        theta, salt,
1535       O        RHO1K,       O        RHO1K,
1536       I        myThid )       I        myThid )
1537    
1538    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1539    CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey
1540    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k, key=kkppkey
1541    CADJ STORE rho1k (:,:)          = comlev1_kpp_k, key=kkppkey
1542    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1543    
1544           call FIND_ALPHA(           call FIND_ALPHA(
1545       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1546       O        WORK1 )       O        WORK1 )
1547    
1548           call FIND_BETA(           call FIND_BETA(
1549       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1550       O        WORK2 )       O        WORK2 )
1551    
1552           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1553              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1554                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1555                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1556                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
# Line 1432  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO Line 1560  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1560              END DO              END DO
1561           END DO           END DO
1562    
1563    #ifdef ALLOW_DIAGNOSTICS
1564             IF ( useDiagnostics ) THEN
1565    c     work1 - density of t(k-1)  & s(k-1) at depth 1
1566    c     work2 - density of t(k  )  & s(k  ) at depth 1
1567    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1568                call FIND_RHO(
1569         I           bi, bj, 1, sNx, 1, sNy, K-1, 1, theta, salt,
1570         O           WORK1,
1571         I           myThid )
1572                call FIND_RHO(
1573         I           bi, bj, 1, sNx, 1, sNy, K  , 1, theta, salt,
1574         O           WORK2,
1575         I           myThid )
1576                DO J = 1, sNy
1577                   DO I = 1, sNx
1578                      IF ( k .LE. klowC(I,J,bi,bj) .AND.
1579         &                 WORK1(I,J) .LT. WORK3(I,J) .AND.
1580         &                 WORK2(I,J) .GE. WORK3(I,J) ) THEN
1581                         KPPMLD(I,J) = MIN ( KPPMLD(I,J),
1582         &                    ABS(((WORK3(I,J)-WORK1(I,J))*rC(k)+
1583         &                    (WORK2(I,J)-WORK3(I,J))*rC(k-1))/
1584         &                    (WORK2(I,J)-WORK1(I,J))))
1585                      ENDIF
1586                   END DO
1587                END DO
1588             ENDIF
1589    #endif /* ALLOW_DIAGNOSTICS */
1590    
1591        END DO        END DO
1592    
1593  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1594        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1595           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1596              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1597              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1598              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1599           END DO           END DO
1600        END DO        END DO
1601    
1602    #ifdef ALLOW_DIAGNOSTICS
1603          IF ( useDiagnostics ) THEN
1604             CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)
1605          ENDIF
1606    #endif /* ALLOW_DIAGNOSTICS */
1607    
1608  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1609    
1610        RETURN        RETURN
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22