/[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.35 by jmc, Fri Oct 19 19:11:17 2007 UTC revision 1.40 by dfer, Thu Oct 23 18:11:00 2008 UTC
# Line 243  C     mixing coefficients blmc to pure R Line 243  C     mixing coefficients blmc to pure R
243                 diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )                 diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
244                 diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )                 diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
245              else              else
246                 ghat(i,k) = 0.                 ghat(i,k) = 0. _d 0
247              endif              endif
248           end do           end do
249        end do        end do
# Line 290  c Line 290  c
290        IMPLICIT NONE        IMPLICIT NONE
291    
292  #include "SIZE.h"  #include "SIZE.h"
293    #include "EEPARAMS.h"
294    #include "PARAMS.h"
295  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
296    
297  c input  c input
# Line 357  c wm, ws    : turbulent velocity scales Line 359  c wm, ws    : turbulent velocity scales
359        _RL         minusone        _RL         minusone
360        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
361    
362    #ifdef ALLOW_DIAGNOSTICS
363    c     KPPBFSFC - Bo+radiation absorbed to d=hbf*hbl + plume (m^2/s^3)
364          _RL KPPBFSFC(imt,Nr)
365          _RL KPPRi(imt,Nr)
366    #endif /* ALLOW_DIAGNOSTICS */
367    
368  c find bulk Richardson number at every grid level until > Ricr  c find bulk Richardson number at every grid level until > Ricr
369  c  c
370  c note: the reference depth is -epsilon/2.*zgrid(k), but the reference  c note: the reference depth is -epsilon/2.*zgrid(k), but the reference
# Line 369  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid Line 377  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid
377  c     initialize hbl and kbl to bottomed out values  c     initialize hbl and kbl to bottomed out values
378    
379        do i = 1, imt        do i = 1, imt
380           Rib(i,1) = 0.0           Rib(i,1) = 0. _d 0
381           kbl(i) = max(kmtj(i),1)           kbl(i) = max(kmtj(i),1)
382           hbl(i) = -zgrid(kbl(i))           hbl(i) = -zgrid(kbl(i))
383        end do        end do
384    
385    #ifdef ALLOW_DIAGNOSTICS
386          do kl = 1, Nr
387             do i = 1, imt
388                KPPBFSFC(i,kl) = 0. _d 0
389                KPPRi(i,kl) = 0. _d 0
390             enddo
391          enddo
392    #endif /* ALLOW_DIAGNOSTICS */
393    
394        do kl = 2, Nr        do kl = 2, Nr
395    
396  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 405  c     compute bfsfc= Bo + radiative cont Line 422  c     compute bfsfc= Bo + radiative cont
422           end do           end do
423  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
424  c     compute bfsfc = plume fraction at hbf * zgrid  c     compute bfsfc = plume fraction at hbf * zgrid
425             IF ( useSALT_PLUME ) THEN
426               do i = 1, imt
427                  worka(i) = zgrid(kl)
428               enddo
429               call SALT_PLUME_FRAC(
430         I         imt, hbf,SPDepth,
431         U         worka,
432         I         myTime, myIter, myThid)
433               do i = 1, imt
434                  bfsfc(i) = bfsfc(i) + boplume(i)*(1. - worka(i))
435               enddo
436             ENDIF
437    #endif /* ALLOW_SALT_PLUME */
438    
439    #ifdef ALLOW_DIAGNOSTICS
440           do i = 1, imt           do i = 1, imt
441              worka(i) = zgrid(kl)              KPPBFSFC(i,kl) = bfsfc(i)
          enddo  
          call PLUMEFRAC(  
      I       imt, hbf,SPDepth,  
      U       worka,  
      I       myTime, myIter, myThid)  
          do i = 1, imt  
             bfsfc(i) = bfsfc(i) + boplume(i)*(1. - worka(i))  
442           enddo           enddo
443  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_DIAGNOSTICS */
444    
445           do i = 1, imt           do i = 1, imt
446              stable(i) = p5 + sign(p5,bfsfc(i))              stable(i) = p5 + sign(p5,bfsfc(i))
# Line 441  c--------------------------------------- Line 466  c---------------------------------------
466       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+
467       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))
468    
469              if (bvsq .eq. 0.) then              if (bvsq .eq. 0. _d 0) then
470                vtsq = 0.0                vtsq = 0. _d 0
471              else              else
472                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc
473              endif              endif
# Line 462  c Line 487  c
487              tempVar1  = dvsq(i,kl) + vtsq                        tempVar1  = dvsq(i,kl) + vtsq          
488              tempVar2 = max(tempVar1, phepsi)              tempVar2 = max(tempVar1, phepsi)
489              Rib(i,kl) = Ritop(i,kl) / tempVar2              Rib(i,kl) = Ritop(i,kl) / tempVar2
490    #ifdef ALLOW_DIAGNOSTICS
491                KPPRi(i,kl) = Rib(i,kl)
492    #endif
493    
494           end do           end do
495        end do        end do
496    
497    #ifdef ALLOW_DIAGNOSTICS
498             CALL DIAGNOSTICS_FILL(KPPBFSFC,'KPPbfsfc',0,Nr,0,1,1,myThid)
499             CALL DIAGNOSTICS_FILL(KPPRi   ,'KPPRi   ',0,Nr,0,1,1,myThid)
500    #endif /* ALLOW_DIAGNOSTICS */
501    
502  cph(  cph(
503  cph  without this store, there is a recomputation error for  cph  without this store, there is a recomputation error for
504  cph  rib in adbldepth (probably partial recomputation problem)      cph  rib in adbldepth (probably partial recomputation problem)    
# Line 516  CADJ &   , key = ikppkey, shape = (/ (sN Line 549  CADJ &   , key = ikppkey, shape = (/ (sN
549        end do        end do
550    
551  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
552        do i = 1, imt        IF ( useSALT_PLUME ) THEN
553           worka(i) = hbl(i)          do i = 1, imt
554        enddo             worka(i) = hbl(i)
555        call PLUMEFRAC(          enddo
556       I       imt,minusone,SPDepth,          call SALT_PLUME_FRAC(
557       U       worka,       I         imt,minusone,SPDepth,
558       I       myTime, myIter, myThid )       U         worka,
559        do i = 1, imt       I         myTime, myIter, myThid )
560           bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))          do i = 1, imt
561        enddo             bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))
562            enddo
563          ENDIF
564  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
565  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
566  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
# Line 601  CADJ &   , key = ikppkey, shape = (/ (sN Line 636  CADJ &   , key = ikppkey, shape = (/ (sN
636        end do        end do
637    
638  #ifdef ALLOW_SALT_PLUME  #ifdef ALLOW_SALT_PLUME
639        do i = 1, imt        IF ( useSALT_PLUME ) THEN
640           worka(i) = hbl(i)          do i = 1, imt
641        enddo             worka(i) = hbl(i)
642        call PLUMEFRAC(          enddo
643       I       imt,minusone,SPDepth,          call SALT_PLUME_FRAC(
644       U       worka,       I         imt,minusone,SPDepth,
645       I       myTime, myIter, myThid )       U         worka,
646        do i = 1, imt       I         myTime, myIter, myThid )
647           bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))          do i = 1, imt
648        enddo             bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))
649            enddo
650          ENDIF
651  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
652  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
653  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
# Line 795  CADJ INIT kpp_ri_tape_mr = common, 1 Line 832  CADJ INIT kpp_ri_tape_mr = common, 1
832  #endif  #endif
833    
834  c constants  c constants
835        c1 = 1.0        c1 = 1. _d 0
836        c0 = 0.0        c0 = 0. _d 0
837    
838  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
839  c     compute interior gradient Ri at all interfaces ki=1,Nr, (not surface)  c     compute interior gradient Ri at all interfaces ki=1,Nr, (not surface)
# Line 843  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 880  CADJ store diffus(:,:,1) = kpp_ri_tape_m
880  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
881    
882           call z121 (           call z121 (
883       U     diffus(1,0,1)       U     diffus(1,0,1),
884       I     myThid )       I     myThid )
885        end do        end do
886  #endif  #endif
# Line 1471  c calculate density, alpha, beta in surf Line 1508  c calculate density, alpha, beta in surf
1508  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1509  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1510  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1511        call FIND_RHO(        CALL FIND_RHO_2D(
1512       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,       I     1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
1513       I     theta, salt,       I     theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1514       O     WORK1,       O     WORK1,
1515       I     myThid )       I     1, bi, bj, myThid )
1516  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1517  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1518  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
# Line 1521  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO Line 1558  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1558  CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1559  CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1560  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1561           call FIND_RHO(           CALL FIND_RHO_2D(
1562       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,       I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1563       I        theta, salt,       I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
1564       O        RHOK,       O        RHOK,
1565       I        myThid )       I        k, bi, bj, myThid )
1566    
1567  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1568  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1569  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1570  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1571           call FIND_RHO(           CALL FIND_RHO_2D(
1572       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K-1, K,       I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1573       I        theta, salt,       I        theta(1-OLx,1-OLy,k-1,bi,bj),salt(1-OLx,1-OLy,k-1,bi,bj),
1574       O        RHOKM1,       O        RHOKM1,
1575       I        myThid )       I        k-1, bi, bj, myThid )
1576    
1577  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1578  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1579  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1580  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1581           call FIND_RHO(           CALL FIND_RHO_2D(
1582       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K,       I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1583       I        theta, salt,       I        theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1584       O        RHO1K,       O        RHO1K,
1585       I        myThid )       I        1, bi, bj, myThid )
1586    
1587  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1588  CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey  CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey
# Line 1577  CADJ STORE rho1k (:,:)          = comlev Line 1614  CADJ STORE rho1k (:,:)          = comlev
1614  c     work1 - density of t(k-1)  & s(k-1) at depth 1  c     work1 - density of t(k-1)  & s(k-1) at depth 1
1615  c     work2 - density of t(k  )  & s(k  ) at depth 1  c     work2 - density of t(k  )  & s(k  ) at depth 1
1616  c     work3 - density of t(1)-.8 & s(1  ) at depth 1  c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1617              call FIND_RHO(              CALL FIND_RHO_2D(
1618       I           bi, bj, 1, sNx, 1, sNy, K-1, 1, theta, salt,       I           1, sNx, 1, sNy, 1,
1619         I           theta(1-OLx,1-OLy,k-1,bi,bj),
1620         I           salt (1-OLx,1-OLy,k-1,bi,bj),
1621       O           WORK1,       O           WORK1,
1622       I           myThid )       I           k-1, bi, bj, myThid )
1623              call FIND_RHO(              CALL FIND_RHO_2D(
1624       I           bi, bj, 1, sNx, 1, sNy, K  , 1, theta, salt,       I           1, sNx, 1, sNy, 1,
1625         I           theta(1-OLx,1-OLy,k,bi,bj),
1626         I           salt (1-OLx,1-OLy,k,bi,bj),
1627       O           WORK2,       O           WORK2,
1628       I           myThid )       I           k, bi, bj, myThid )
1629              DO J = 1, sNy              DO J = 1, sNy
1630                 DO I = 1, sNx                 DO I = 1, sNx
1631                    IF ( k .LE. klowC(I,J,bi,bj) .AND.                    IF ( k .LE. klowC(I,J,bi,bj) .AND.
# Line 1614  c     compute arrays for K = Nrp1 Line 1655  c     compute arrays for K = Nrp1
1655  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
1656        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
1657           CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)           CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)
1658             CALL DIAGNOSTICS_FILL(DBSFC ,'KPPdbsfc',0,Nr,0,1,1,myThid)
1659             CALL DIAGNOSTICS_FILL(DBLOC ,'KPPdbloc',0,Nr,0,1,1,myThid)
1660        ENDIF        ENDIF
1661  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
1662    
Legend:
Removed from v.1.35  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.22