--- trunk/libf/dyn3d/dissip.f 2008/04/18 14:45:53 10 +++ trunk/libf/dyn3d/Dissipation/dissip.f90 2012/01/30 12:54:02 57 @@ -1,143 +1,90 @@ -! -! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/dissip.F,v 1.1.1.1 2004/05/19 12:53:05 lmdzadmin Exp $ -! - SUBROUTINE dissip( vcov,ucov,teta,p, dv,du,dh ) -c - use dimens_m - use paramet_m - use comconst - use comdissnew - use comgeom - IMPLICIT NONE - - -c .. Avec nouveaux operateurs star : gradiv2 , divgrad2, nxgraro2 ... -c ( 10/01/98 ) - -c======================================================================= -c -c Auteur: P. Le Van -c ------- -c -c Objet: -c ------ -c -c Dissipation horizontale -c -c======================================================================= -c----------------------------------------------------------------------- -c Declarations: -c ------------- - - include "comdissipn.h" - -c Arguments: -c ---------- - - REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm) - REAL, intent(in):: p( ip1jmp1,llmp1 ) - REAL dv(ip1jm,llm),du(ip1jmp1,llm),dh(ip1jmp1,llm) - -c Local: -c ------ - - REAL gdx(ip1jmp1,llm),gdy(ip1jm,llm) - REAL grx(ip1jmp1,llm),gry(ip1jm,llm) - REAL te1dt(llm),te2dt(llm),te3dt(llm) - REAL deltapres(ip1jmp1,llm) - - INTEGER l,ij - - REAL SSUM - -c----------------------------------------------------------------------- -c initialisations: -c ---------------- - - DO l=1,llm - te1dt(l) = tetaudiv(l) * dtdiss - te2dt(l) = tetaurot(l) * dtdiss - te3dt(l) = tetah(l) * dtdiss - ENDDO - du=0. - dv=0. - dh=0. - -c----------------------------------------------------------------------- -c Calcul de la dissipation: -c ------------------------- - -c Calcul de la partie grad ( div ) : -c ------------------------------------- - - - IF(lstardis) THEN - CALL gradiv2( llm,ucov,vcov,nitergdiv,gdx,gdy ) - ELSE - CALL gradiv ( llm,ucov,vcov,nitergdiv,gdx,gdy ) - ENDIF - - DO l=1,llm - - DO ij = 1, iip1 - gdx( ij ,l) = 0. - gdx(ij+ip1jm,l) = 0. - ENDDO - - DO ij = iip2,ip1jm - du(ij,l) = du(ij,l) - te1dt(l) *gdx(ij,l) - ENDDO - DO ij = 1,ip1jm - dv(ij,l) = dv(ij,l) - te1dt(l) *gdy(ij,l) - ENDDO - - ENDDO - -c calcul de la partie n X grad ( rot ): -c --------------------------------------- - - IF(lstardis) THEN - CALL nxgraro2( llm,ucov, vcov, nitergrot,grx,gry ) - ELSE - CALL nxgrarot( llm,ucov, vcov, nitergrot,grx,gry ) - ENDIF - - - DO l=1,llm - DO ij = 1, iip1 - grx(ij,l) = 0. - ENDDO - - DO ij = iip2,ip1jm - du(ij,l) = du(ij,l) - te2dt(l) * grx(ij,l) - ENDDO - DO ij = 1, ip1jm - dv(ij,l) = dv(ij,l) - te2dt(l) * gry(ij,l) - ENDDO - ENDDO - -c calcul de la partie div ( grad ): -c ----------------------------------- - - - IF(lstardis) THEN - - DO l = 1, llm - DO ij = 1, ip1jmp1 - deltapres(ij,l) = AMAX1( 0., p(ij,l) - p(ij,l+1) ) - ENDDO - ENDDO - - CALL divgrad2( llm,teta, deltapres ,niterh, gdx ) - ELSE - CALL divgrad ( llm,teta, niterh, gdx ) - ENDIF - - DO l = 1,llm - DO ij = 1,ip1jmp1 - dh( ij,l ) = dh( ij,l ) - te3dt(l) * gdx( ij,l ) - ENDDO - ENDDO +module dissip_m - RETURN - END + IMPLICIT NONE + +contains + + SUBROUTINE dissip(vcov, ucov, teta, p, dv, du, dh) + + ! From dyn3d/dissip.F, version 1.1.1.1 2004/05/19 12:53:05 + ! Author: P. Le Van + ! Objet : calcul de la dissipation horizontale + ! Avec opérateurs star : gradiv2, divgrad2, nxgraro2 + + USE dimens_m, ONLY: iim, jjm, llm + USE comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh + USE inidissip_m, ONLY: dtdiss, tetah, tetaudiv, tetaurot, cdivu, crot, cdivh + use gradiv2_m, only: gradiv2 + use nr_util, only: assert + + REAL, intent(in):: vcov(:, :, :) ! (iim + 1, jjm, llm) + REAL, intent(in):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) + REAL, intent(in):: teta(:, :, :) ! (iim + 1, jjm + 1, llm) + REAL, INTENT(IN):: p(:, :, :) ! (iim + 1, jjm + 1, llm + 1) + REAL, intent(out):: dv(:, :, :) ! (iim + 1, jjm, llm) + REAL, intent(out):: du(:, :, :) ! (iim + 1, jjm + 1, llm) + REAL, intent(out):: dh(:, :, :) ! (iim + 1, jjm + 1, llm) + + ! Local: + REAL gdx(iim + 1, jjm + 1, llm), gdy(iim + 1, jjm, llm) + REAL grx(iim + 1, jjm + 1, llm), gry(iim + 1, jjm, llm) + REAL tedt(llm) + REAL deltapres(iim + 1, jjm + 1, llm) + INTEGER l + + !----------------------------------------------------------------------- + + call assert((/size(vcov, 1), size(ucov, 1), size(teta, 1), size(p, 1), & + size(dv, 1), size(du, 1), size(dh, 1)/) == iim + 1, "dissip iim") + call assert((/size(vcov, 2), size(ucov, 2) - 1, size(teta, 2) - 1, & + size(p, 2) - 1, size(dv, 2), size(du, 2) - 1, size(dh, 2) - 1/) & + == jjm, "dissip jjm") + call assert((/size(vcov, 3), size(ucov, 3), size(teta, 3), size(p, 3) - 1, & + size(dv, 3), size(du, 3), size(dh, 3)/) == llm, "dissip llm") + + du(:, 1, :) = 0. + du(:, jjm + 1, :) = 0. + + ! Calcul de la partie grad (div) : + + IF (lstardis) THEN + CALL gradiv2(ucov, vcov, nitergdiv, gdx, gdy, cdivu) + ELSE + CALL gradiv(llm, ucov, vcov, nitergdiv, gdx, gdy, cdivu) + END IF + + tedt = tetaudiv * dtdiss + forall (l = 1: llm) + du(:, 2: jjm, l) = - tedt(l) * gdx(:, 2: jjm, l) + dv(:, :, l) = - tedt(l) * gdy(:, :, l) + END forall + + ! Calcul de la partie n X grad (rot) : + + IF (lstardis) THEN + CALL nxgraro2(llm, ucov, vcov, nitergrot, grx, gry, crot) + ELSE + CALL nxgrarot(llm, ucov, vcov, nitergrot, grx, gry, crot) + END IF + + tedt = tetaurot * dtdiss + forall (l = 1: llm) + du(:, 2: jjm, l) = du(:, 2: jjm, l) - tedt(l) * grx(:, 2: jjm, l) + dv(:, :, l) = dv(:, :, l) - tedt(l) * gry(:, :, l) + END forall + + ! calcul de la partie div (grad) : + + IF (lstardis) THEN + forall (l = 1: llm) & + deltapres(:, :, l) = max(0., p(:, :, l) - p(:, :, l + 1)) + CALL divgrad2(llm, teta, deltapres, niterh, gdx, cdivh) + ELSE + CALL divgrad(llm, teta, niterh, gdx, cdivh) + END IF + + forall (l = 1: llm) dh(:, :, l) = - tetah(l) * dtdiss * gdx(:, :, l) + + END SUBROUTINE dissip + +end module dissip_m