/[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.1 by adcroft, Wed Jun 21 19:45:52 2000 UTC revision 1.2 by heimbach, Tue Sep 12 18:14:32 2000 UTC
# Line 152  c--------------------------------------- Line 152  c---------------------------------------
152       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
153       &     )       &     )
154    
155  CADJ STORE hbl,bfsfc,stable,casea,kbl = kpptape, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey
156    
157  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
158  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 163  c--------------------------------------- Line 163  c---------------------------------------
163       O     , dkm1, blmc, ghat, sigma       O     , dkm1, blmc, ghat, sigma
164       &     )       &     )
165    
166  CADJ STORE dkm1,blmc,ghat = kpptape, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey
167    
168  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
169  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 281  c zwork     : depth work array Line 281  c zwork     : depth work array
281        _RS wm(imt), ws(imt)        _RS wm(imt), ws(imt)
282        _RS zwork(imt)        _RS zwork(imt)
283    
284        _RS bvsq, vtsq, hekman, hmonob, hlimit, p5        _RS bvsq, vtsq, hekman, hmonob, hlimit
285        integer i, kl        integer i, kl
286    
287        parameter (p5 = 0.5)        _RS        p5    , eins    , m1
288          parameter (p5=0.5, eins=1.0, m1=-1.0 )
289    
290  crg intermediate result in RL to avoid overflow in adjoint  crg intermediate result in RL to avoid overflow in adjoint
291         _RL        dpshear2         _RL        dpshear2
# Line 388  ctl end-replace Line 389  ctl end-replace
389           end do           end do
390        end do        end do
391    
392  CADJ store kbl = kpptape  CADJ store kbl = comlev1_kpp
393  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
394    
395        do i = 1, imt        do i = 1, imt
# Line 408  c     linearly interpolate to find hbl w Line 409  c     linearly interpolate to find hbl w
409           endif           endif
410        end do        end do
411    
412  CADJ store hbl = kpptape  CADJ store hbl = comlev1_kpp
413  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
414    
415  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 416  c     find stability and buoyancy forcin Line 417  c     find stability and buoyancy forcin
417  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
418    
419        call SWFRAC(        call SWFRAC(
420       I     imt, -1.0, hbl,       I     imt, m1, hbl,
421       O     bfsfc)       O     bfsfc)
422    
423    CADJ store bfsfc = comlev1_kpp
424    CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
425    
426        do i = 1, imt        do i = 1, imt
427           bfsfc(i)  = bo(i) + bosol(i) * (1. - bfsfc(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - bfsfc(i))
428           stable(i) = p5 + sign( p5, bfsfc(i) )           stable(i) = p5 + sign( p5, bfsfc(i) )
429    
430  #ifdef KPP_TEST_DENOM  #ifdef KPP_TEST_DENOM
431  ctl add  ctl add
432          bfsfc(i) = sign(1.0,bfsfc(i))*max(phepsi,abs(bfsfc(i)))          bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
433  ctl end-add  ctl end-add
434  #else  #else
435  c--      ensures bfsfc is never 0  c--      ensures bfsfc is never 0
# Line 458  ctl end-replace Line 462  ctl end-replace
462           kbl(i) = kmtj(i)           kbl(i) = kmtj(i)
463        end do        end do
464    
465  CADJ store hbl = kpptape  CADJ store hbl = comlev1_kpp
466  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
467    
468  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 478  c      find stability and buoyancy forci Line 482  c      find stability and buoyancy forci
482  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
483    
484        call SWFRAC(        call SWFRAC(
485       I     imt, -1.0, hbl,       I     imt, m1, hbl,
486       O     bfsfc)       O     bfsfc)
487    
488    CADJ store bfsfc = comlev1_kpp
489    CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
490    
491        do i = 1, imt        do i = 1, imt
492           bfsfc(i) = bo(i) + bosol(i) * (1. - bfsfc(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - bfsfc(i))
493           stable(i) = p5 + sign( p5, bfsfc(i) )           stable(i) = p5 + sign( p5, bfsfc(i) )
494  #ifdef KPP_TEST_DENOM  #ifdef KPP_TEST_DENOM
495  ctl add  ctl add
496          bfsfc(i) = sign(1.0,bfsfc(i))*max(phepsi,abs(bfsfc(i)))          bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
497  ctl end-add  ctl end-add
498  #else  #else
499  c--      ensures bfsfc is never 0  c--      ensures bfsfc is never 0
# Line 685  c--------------------------------------- Line 692  c---------------------------------------
692  c     vertically smooth Ri  c     vertically smooth Ri
693    
694        do mr = 1, num_v_smooth_Ri        do mr = 1, num_v_smooth_Ri
695    
696    CADJ store diffus(:,:,1) = comlev1_kpp_sm
697    CADJ &  , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
698    
699           call z121 (           call z121 (
700       U     diffus(1,0,1))       U     diffus(1,0,1))
701        end do        end do
702    
703  CADJ store diffus = kpptape  CADJ store diffus = comlev1_kpp
704  CADJ &  , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2,3 /)  CADJ &  , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2,3 /)
705    
706  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 718  c    mixing due to internal waves, and s Line 729  c    mixing due to internal waves, and s
729    
730              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
731              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0
732              diffus(i,ki,3) = diffus(i,ki,2)              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0
733    
734           end do           end do
735        end do        end do
# Line 746  c     Apply 121 smoothing in k to 2-d ar Line 757  c     Apply 121 smoothing in k to 2-d ar
757  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
758  c     bottom (Nrp1) value is set to input value from above.  c     bottom (Nrp1) value is set to input value from above.
759    
760  c     Note that its important to exclude from the smoothing any points  c     Note that it is important to exclude from the smoothing any points
761  c     that are outside the range of the K(Ri) scheme, ie.  >0.8, or <0.0.  c     that are outside the range of the K(Ri) scheme, ie.  >0.8, or <0.0.
762  c     Otherwise, there is interfence with other physics, especially  c     Otherwise, there is interference with other physics, especially
763  c     penetrative convection.  c     penetrative convection.
764    
765        IMPLICIT NONE        IMPLICIT NONE
# Line 764  c v     : 2-D array to be smoothed in Nr Line 775  c v     : 2-D array to be smoothed in Nr
775    
776  c local  c local
777        _RS zwork, zflag        _RS zwork, zflag
778        _RS KRi_range(0:Nrp1)        _RS KRi_range(1:Nrp1)
779        integer i, k, km1, kp1        integer i, k, km1, kp1
780    
781        _RS         p25       , p5      , p2        _RS         p0      , p25       , p5      , p2
782        parameter ( p25 = 0.25, p5 = 0.5, p2 = 2.0 )        parameter ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 )
783    
784          KRi_range(Nrp1) = p0
785    
786    #ifdef ALLOW_AUTODIFF_TAMC
787    C--   dummy assignment to end declaration part for TAMC
788          i = 0
789    
790    C--   HPF directive to help TAMC
791    CHPF$ INDEPENDENT
792    #endif /* ALLOW_AUTODIFF_TAMC */
793        do i = 1, imt        do i = 1, imt
794    
795           v(i,Nrp1) = v(i,Nr)           v(i,Nrp1) = v(i,Nr)
# Line 786  c local Line 806  c local
806           zflag  = p2 + KRi_range(1) * KRi_range(2)           zflag  = p2 + KRi_range(1) * KRi_range(2)
807           v(i,1) = v(i,1) / zflag           v(i,1) = v(i,1) / zflag
808    
809    CADJ INIT z121tape = common, Nr
810           do k = 2, Nr           do k = 2, Nr
811    CADJ STORE v(i,k), zwork = z121tape
812              km1 = k - 1              km1 = k - 1
813              kp1 = k + 1              kp1 = k + 1
814              zflag = v(i,k)              zflag = v(i,k)
# Line 1060  c     sigma(imt)                  normal Line 1082  c     sigma(imt)                  normal
1082  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1083    
1084  c  local  c  local
1085  c     gat1  (imt,mdiff)          shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1086  c     dat1  (imt,mdiff)          derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1087  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1088        _RS gat1  (imt,mdiff)        _RS gat1m(imt), gat1s(imt), gat1t(imt)
1089        _RS dat1  (imt,mdiff)        _RS dat1m(imt), dat1s(imt), dat1t(imt)
1090        _RS ws(imt), wm(imt)        _RS ws(imt), wm(imt)
1091        integer i, kn, ki        integer i, kn, ki
1092        _RS     R, dvdzup, dvdzdn, viscp        _RS     R, dvdzup, dvdzdn, viscp
1093        _RS     difsp, diftp, visch, difsh, difth        _RS     difsp, diftp, visch, difsh, difth
1094        _RS     f1, sig, a1, a2, a3, delhat        _RS     f1, sig, a1, a2, a3, delhat
1095        _RS     GM, Gs, Gt, p0        _RS     Gm, Gs, Gt
1096          _RL     dum
1097    
1098        parameter (p0 = 0.0)        _RS        p0    , eins
1099          parameter (p0=0.0, eins=1.0)
       _RL dum  
1100    
1101  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1102  c compute velocity scales at hbl  c compute velocity scales at hbl
# Line 1122  c--------------------------------------- Line 1144  c---------------------------------------
1144  ctl replace (Important!!! not phepsi**4 !!!)  ctl replace (Important!!! not phepsi**4 !!!)
1145           f1 = stable(i) * conc1 * bfsfc(i) /           f1 = stable(i) * conc1 * bfsfc(i) /
1146       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1147           wm(i) = sign(1.0,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1148           gat1(i,1) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1149           dat1(i,1) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
1150  #else  #else
1151           f1 = stable(i) * conc1 * bfsfc(i) / (ustar(i)**4+epsln)           f1 = stable(i) * conc1 * bfsfc(i) / (ustar(i)**4+epsln)
1152           gat1(i,1) = visch / hbl(i) / (wm(i)+epsln)           gat1m(i) = visch / hbl(i) / (wm(i)+epsln)
1153           dat1(i,1) = -viscp / (wm(i)+epsln) + f1 * visch           dat1m(i) = -viscp / (wm(i)+epsln) + f1 * visch
1154  #endif  #endif
1155           dat1(i,1) = min(dat1(i,1),p0)           dat1m(i) = min(dat1m(i),p0)
1156    
1157  #ifdef KPP_TEST_DENOM  #ifdef KPP_TEST_DENOM
1158           ws(i) = sign(1.0,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1159           gat1(i,2) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1160           dat1(i,2) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
1161  #else  #else
1162           gat1(i,2) = difsh / hbl(i) / (ws(i)+epsln)           gat1s(i) = difsh / hbl(i) / (ws(i)+epsln)
1163           dat1(i,2) = -difsp / (ws(i)+epsln) + f1 * difsh           dat1s(i) = -difsp / (ws(i)+epsln) + f1 * difsh
1164  #endif  #endif
1165           dat1(i,2) = min(dat1(i,2),p0)           dat1s(i) = min(dat1s(i),p0)
1166    
1167  #ifdef KPP_TEST_DENOM  #ifdef KPP_TEST_DENOM
1168           gat1(i,3) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1169           dat1(i,3) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
1170  #else  #else
1171           gat1(i,3) = difth / hbl(i) / (ws(i)+epsln)           gat1t(i) = difth / hbl(i) / (ws(i)+epsln)
1172           dat1(i,3) = -diftp / (ws(i)+epsln) + f1 * difth           dat1t(i) = -diftp / (ws(i)+epsln) + f1 * difth
1173  #endif  #endif
1174           dat1(i,3) = min(dat1(i,3),p0)           dat1t(i) = min(dat1t(i),p0)
1175    
1176        end do        end do
1177    
# Line 1177  c--------------------------------------- Line 1199  c---------------------------------------
1199              a2 = 3. - 2. * sig              a2 = 3. - 2. * sig
1200              a3 = sig - 1.              a3 = sig - 1.
1201    
1202              GM = a1 + a2 * gat1(i,1) + a3 * dat1(i,1)              Gm = a1 + a2 * gat1m(i) + a3 * dat1m(i)
1203              Gs = a1 + a2 * gat1(i,2) + a3 * dat1(i,2)              Gs = a1 + a2 * gat1s(i) + a3 * dat1s(i)
1204              Gt = a1 + a2 * gat1(i,3) + a3 * dat1(i,3)              Gt = a1 + a2 * gat1t(i) + a3 * dat1t(i)
1205    
1206  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1207  c     compute boundary layer diffusivities at the interfaces  c     compute boundary layer diffusivities at the interfaces
1208  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1209    
1210              blmc(i,ki,1) = hbl(i) * wm(i) * sig * (1. + sig * GM)              blmc(i,ki,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)
1211              blmc(i,ki,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)              blmc(i,ki,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)
1212              blmc(i,ki,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)              blmc(i,ki,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)
1213    
# Line 1224  c--------------------------------------- Line 1246  c---------------------------------------
1246           a1 = sig - 2.           a1 = sig - 2.
1247           a2 = 3. - 2. * sig           a2 = 3. - 2. * sig
1248           a3 = sig - 1.           a3 = sig - 1.
1249           GM = a1 + a2 * gat1(i,1) + a3 * dat1(i,1)           Gm = a1 + a2 * gat1m(i) + a3 * dat1m(i)
1250           Gs = a1 + a2 * gat1(i,2) + a3 * dat1(i,2)           Gs = a1 + a2 * gat1s(i) + a3 * dat1s(i)
1251           Gt = a1 + a2 * gat1(i,3) + a3 * dat1(i,3)           Gt = a1 + a2 * gat1t(i) + a3 * dat1t(i)
1252           dkm1(i,1) = hbl(i) * wm(i) * sig * (1. + sig * GM)           dkm1(i,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)
1253           dkm1(i,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)           dkm1(i,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)
1254           dkm1(i,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)           dkm1(i,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)
1255        end do        end do
# Line 1417  c calculate density, alpha, beta in surf Line 1439  c calculate density, alpha, beta in surf
1439    
1440  c calculate alpha, beta, and gradients in interior layers  c calculate alpha, beta, and gradients in interior layers
1441    
1442  !HPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1443        DO K = 2, Nr        DO K = 2, Nr
1444    
1445           call FIND_RHO(           call FIND_RHO(
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22