--- trunk/libf/dyn3d/inidissip.f 2010/03/05 16:43:45 25 +++ trunk/libf/dyn3d/inidissip.f90 2010/03/09 15:27:15 26 @@ -1,218 +1,197 @@ -! -! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/inidissip.F,v 1.1.1.1 2004/05/19 12:53:06 lmdzadmin Exp $ -! - SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , - * tetagdiv,tetagrot,tetatemp ) -c======================================================================= -c initialisation de la dissipation horizontale -c======================================================================= -c----------------------------------------------------------------------- -c declarations: -c ------------- - - use dimens_m - use paramet_m - use comconst - use comvert - use conf_gcm_m - use comdissipn - IMPLICIT NONE - - LOGICAL lstardis - INTEGER nitergdiv,nitergrot,niterh - REAL tetagdiv,tetagrot,tetatemp - REAL fact,zvert(llm),zz - REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm) - REAL ullm,vllm,umin,vmin,zhmin,zhmax - REAL zllm,z1llm - - INTEGER l,ij,idum,ii - REAL tetamin - - REAL ran1 - - -c----------------------------------------------------------------------- - - print *, "Call sequence information: inidissip" -c -c calcul des valeurs propres des operateurs par methode iterrative: -c ----------------------------------------------------------------- - - crot = -1. - cdivu = -1. - cdivh = -1. - -c calcul de la valeur propre de divgrad: -c -------------------------------------- - idum = 0 - DO l = 1, llm +module inidissip_m + + use dimens_m, only: llm + + IMPLICIT NONE + + private llm + + REAL dtdiss + integer idissip ! periode de la dissipation (en pas) + real tetaudiv(llm),tetaurot(llm),tetah(llm) + real cdivu, crot, cdivh + +contains + + SUBROUTINE inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, & + tetagrot, tetatemp) + + ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06 + ! Initialisation de la dissipation horizontale + + USE comconst, ONLY : dtvr + USE comvert, ONLY : preff, presnivs + USE conf_gcm_m, ONLY : iperiod + USE dimens_m, ONLY : jjm, llm + USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1 + use new_unit_m, only: new_unit + + LOGICAL, intent(in):: lstardis + INTEGER, intent(in):: nitergdiv, nitergrot, niterh + REAL, intent(in):: tetagdiv, tetagrot, tetatemp + + ! Variables local to the procedure: + REAL zvert(llm) + REAL zh(ip1jmp1), zu(ip1jmp1), zv(ip1jm), deltap(ip1jmp1, llm) + REAL ullm, vllm, umin, vmin, zhmin, zhmax + REAL zllm, z1llm + INTEGER l, ij, idum, ii, unit + REAL tetamin + REAL ran1 + + !----------------------------------------------------------------------- + + PRINT *, 'Call sequence information: inidissip' + + ! calcul des valeurs propres des operateurs par methode iterrative: + + crot = -1. + cdivu = -1. + cdivh = -1. + + ! calcul de la valeur propre de divgrad: + + idum = 0 + DO l = 1, llm DO ij = 1, ip1jmp1 - deltap(ij,l) = 1. - ENDDO - ENDDO - - idum = -1 - zh(1) = RAN1(idum)-.5 - idum = 0 - DO ij = 2, ip1jmp1 - zh(ij) = RAN1(idum) -.5 - ENDDO - - CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) - - CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) - - IF ( zhmin .GE. zhmax ) THEN - PRINT*,' Inidissip zh min max ',zhmin,zhmax - STOP'probleme generateur alleatoire dans inidissip' - ENDIF - - zllm = ABS( zhmax ) - DO l = 1,50 - IF(lstardis) THEN - CALL divgrad2(1,zh,deltap,niterh,zh) - ELSE - CALL divgrad (1,zh,niterh,zh) - ENDIF - - CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) - - zllm = ABS( zhmax ) - z1llm = 1./zllm - DO ij = 1,ip1jmp1 - zh(ij) = zh(ij)* z1llm - ENDDO - ENDDO - - IF(lstardis) THEN - cdivh = 1./ zllm - ELSE - cdivh = zllm ** ( -1./niterh ) - ENDIF - -c calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) -c ----------------------------------------------------------------- - print*,'calcul des valeurs propres' - - DO 20 ii = 1, 2 -c - DO ij = 1, ip1jmp1 - zu(ij) = RAN1(idum) -.5 - ENDDO - CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) - DO ij = 1, ip1jm - zv(ij) = RAN1(idum) -.5 - ENDDO - CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) - - CALL minmax(iip1*jjp1,zu,umin,ullm ) - CALL minmax(iip1*jjm, zv,vmin,vllm ) - - ullm = ABS ( ullm ) - vllm = ABS ( vllm ) - - DO 5 l = 1, 50 - IF(ii.EQ.1) THEN -ccccc CALL covcont( 1,zu,zv,zu,zv ) - IF(lstardis) THEN - CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv ) - ELSE - CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv ) - ENDIF - ELSE - IF(lstardis) THEN - CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv ) - ELSE - CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv ) - ENDIF - ENDIF - - CALL minmax(iip1*jjp1,zu,umin,ullm ) - CALL minmax(iip1*jjm, zv,vmin,vllm ) - - ullm = ABS ( ullm ) - vllm = ABS ( vllm ) - - zllm = MAX( ullm,vllm ) - z1llm = 1./ zllm - DO ij = 1, ip1jmp1 - zu(ij) = zu(ij)* z1llm - ENDDO - DO ij = 1, ip1jm - zv(ij) = zv(ij)* z1llm - ENDDO - 5 CONTINUE - - IF ( ii.EQ.1 ) THEN - IF(lstardis) THEN - cdivu = 1./zllm - ELSE - cdivu = zllm **( -1./nitergdiv ) - ENDIF - ELSE - IF(lstardis) THEN - crot = 1./ zllm - ELSE - crot = zllm **( -1./nitergrot ) - ENDIF - ENDIF - - 20 CONTINUE - -c petit test pour les operateurs non star: -c ---------------------------------------- - -c IF(.NOT.lstardis) THEN - fact = rad*24./float(jjm) - fact = fact*fact - PRINT*,'coef u ', fact/cdivu, 1./cdivu - PRINT*,'coef r ', fact/crot , 1./crot - PRINT*,'coef h ', fact/cdivh, 1./cdivh -c ENDIF - -c----------------------------------------------------------------------- -c variation verticale du coefficient de dissipation: -c -------------------------------------------------- - - DO l=1,llm - zvert(l)=1. - ENDDO - - fact=2. -c - DO l = 1, llm - zz = 1. - preff/presnivs(l) - zvert(l)= fact -( fact-1.)/( 1.+zz*zz ) - ENDDO - - - PRINT*,'Constantes de temps de la diffusion horizontale' - - tetamin = 1.e+6 - - DO l=1,llm - tetaudiv(l) = zvert(l)/tetagdiv - tetaurot(l) = zvert(l)/tetagrot - tetah(l) = zvert(l)/tetatemp - - IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) - IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) - IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) - ENDDO - - PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod - idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod - PRINT *,' tetamin = ',tetamin - idissip = MAX(iperiod,idissip) - dtdiss = idissip * dtvr - PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss - - DO l = 1,llm - PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), - * dtdiss*tetah(l) - ENDDO - -c - RETURN - END + deltap(ij, l) = 1. + END DO + END DO + + idum = -1 + zh(1) = ran1(idum) - .5 + idum = 0 + DO ij = 2, ip1jmp1 + zh(ij) = ran1(idum) - .5 + END DO + + CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1) + + CALL minmax(iip1*jjp1, zh, zhmin, zhmax) + + IF (zhmin>=zhmax) THEN + PRINT *, ' Inidissip zh min max ', zhmin, zhmax + STOP 'probleme generateur alleatoire dans inidissip' + END IF + + zllm = abs(zhmax) + DO l = 1, 50 + IF (lstardis) THEN + CALL divgrad2(1, zh, deltap, niterh, zh) + ELSE + CALL divgrad(1, zh, niterh, zh) + END IF + + CALL minmax(iip1*jjp1, zh, zhmin, zhmax) + + zllm = abs(zhmax) + z1llm = 1./zllm + DO ij = 1, ip1jmp1 + zh(ij) = zh(ij)*z1llm + END DO + END DO + + IF (lstardis) THEN + cdivh = 1./zllm + ELSE + cdivh = zllm**(-1./niterh) + END IF + + ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) + + PRINT *, 'calcul des valeurs propres' + + DO ii = 1, 2 + + DO ij = 1, ip1jmp1 + zu(ij) = ran1(idum) - .5 + END DO + CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE., 1) + DO ij = 1, ip1jm + zv(ij) = ran1(idum) - .5 + END DO + CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1) + + CALL minmax(iip1*jjp1, zu, umin, ullm) + CALL minmax(iip1*jjm, zv, vmin, vllm) + + ullm = abs(ullm) + vllm = abs(vllm) + + DO l = 1, 50 + IF (ii==1) THEN + IF (lstardis) THEN + CALL gradiv2(1, zu, zv, nitergdiv, zu, zv) + ELSE + CALL gradiv(1, zu, zv, nitergdiv, zu, zv) + END IF + ELSE + IF (lstardis) THEN + CALL nxgraro2(1, zu, zv, nitergrot, zu, zv) + ELSE + CALL nxgrarot(1, zu, zv, nitergrot, zu, zv) + END IF + END IF + + CALL minmax(iip1*jjp1, zu, umin, ullm) + CALL minmax(iip1*jjm, zv, vmin, vllm) + + ullm = abs(ullm) + vllm = abs(vllm) + + zllm = max(ullm, vllm) + z1llm = 1./zllm + DO ij = 1, ip1jmp1 + zu(ij) = zu(ij)*z1llm + END DO + DO ij = 1, ip1jm + zv(ij) = zv(ij)*z1llm + END DO + end DO + + IF (ii==1) THEN + IF (lstardis) THEN + cdivu = 1./zllm + ELSE + cdivu = zllm**(-1./nitergdiv) + END IF + ELSE + IF (lstardis) THEN + crot = 1./zllm + ELSE + crot = zllm**(-1./nitergrot) + END IF + END IF + + END DO + + PRINT *, 'cdivu = ', cdivu + PRINT *, 'crot = ', crot + PRINT *, 'cdivh = ', cdivh + + ! Variation verticale du coefficient de dissipation : + zvert = 2. - 1. / (1. + (1. - preff / presnivs)**2) + + tetaudiv = zvert / tetagdiv + tetaurot = zvert / tetagrot + tetah = zvert / tetatemp + call new_unit(unit) + open(unit, file="inidissip.csv", status="replace", action="write") + write(unit, fmt=*) "tetaudiv tetaurot tetah" ! title line + do l = 1, llm + write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l) + end do + close(unit) + print *, 'Created file "inidissip.csv".' + + tetamin = min(1E6, minval(1. / tetaudiv), minval(1. / tetaurot), & + minval(1. / tetah)) + PRINT *, 'tetamin = ', tetamin + idissip = max(iperiod, int(tetamin / (2 * dtvr * iperiod)) * iperiod) + PRINT *, 'idissip = ', idissip + dtdiss = idissip * dtvr + PRINT *, 'dtdiss = ', dtdiss + + END SUBROUTINE inidissip + +end module inidissip_m