--- trunk/libf/dyn3d/inidissip.f 2010/03/05 16:43:45 25 +++ trunk/dyn3d/Dissipation/inidissip.f90 2013/11/15 18:45:49 76 @@ -1,218 +1,134 @@ -! -! $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 - 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 +module inidissip_m + + use dimens_m, only: llm + + IMPLICIT NONE + + private llm + + REAL dtdiss ! in s + integer idissip ! période de la dissipation (en pas de temps) + real tetaudiv(llm), tetaurot(llm), tetah(llm) ! in s-1 + real cdivu, crot, cdivh + +contains + + SUBROUTINE inidissip + + ! From dyn3d/inidissip.F, version 1.1.1.1 2004/05/19 12:53:06 + + ! Initialisation de la dissipation horizontale. Calcul des valeurs + ! propres des opérateurs par méthode itérative. + + USE comconst, ONLY: dtvr + use comdissnew, only: nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, & + tetatemp + USE disvert_m, ONLY: preff, presnivs + USE conf_gcm_m, ONLY: iperiod + USE dimens_m, ONLY: iim, jjm + use divgrad2_m, only: divgrad2 + use filtreg_m, only: filtreg + use gradiv2_m, only: gradiv2 + use jumble, only: new_unit + use nxgraro2_m, only: nxgraro2 + USE paramet_m, ONLY: jjp1 + + ! Variables local to the procedure: + REAL zvert(llm), max_zvert ! no dimension + REAL, dimension(iim + 1, jjm + 1, 1):: zh, zu, gx, divgra, deltap + real zv(iim + 1, jjm, 1), gy(iim + 1, jjm, 1) + REAL zllm + INTEGER l, seed_size, ii, unit + REAL tetamin ! in s + + !----------------------------------------------------------------------- + + PRINT *, 'Call sequence information: inidissip' + call random_seed(size=seed_size) + call random_seed(put=(/(1, ii = 1, seed_size)/)) + + PRINT *, 'Calcul des valeurs propres de divgrad' + deltap = 1. + call random_number(zh) + zh = zh - 0.5 + CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE.) + + DO l = 1, 50 + CALL divgrad2(1, zh, deltap, niterh, divgra, -1.) + zllm = abs(maxval(divgra)) + zh = divgra / zllm + END DO + + cdivh = 1. / zllm + PRINT *, 'cdivh = ', cdivh + + PRINT *, 'Calcul des valeurs propres de gradiv' + call random_number(zu) + zu = zu - 0.5 + CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.) + call random_number(zv) + zv = zv - 0.5 + CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.) + + DO l = 1, 50 + CALL gradiv2(zu, zv, nitergdiv, gx, gy, -1.) + zllm = max(abs(maxval(gx)), abs(maxval(gy))) + zu = gx / zllm + zv = gy / zllm + end DO + + cdivu = 1. / zllm + PRINT *, 'cdivu = ', cdivu + + PRINT *, 'Calcul des valeurs propres de nxgrarot' + call random_number(zu) + zu = zu - 0.5 + CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.) + call random_number(zv) + zv = zv - 0.5 + CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.) + + DO l = 1, 50 + CALL nxgraro2(zu, zv, nitergrot, gx, gy, -1.) + zllm = max(abs(maxval(gx)), abs(maxval(gy))) + zu = gx / zllm + zv = gy / zllm + end DO + + crot = 1. / zllm + PRINT *, 'crot = ', crot + + ! Variation verticale du coefficient de dissipation : + zvert = 2. - 1. / (1. + (preff / presnivs - 1.)**2) + ! (between 1 and 2) + + tetaudiv = zvert / tetagdiv + tetaurot = zvert / tetagrot + tetah = zvert / tetatemp + + max_zvert = maxval(zvert) + tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, & + tetatemp / max_zvert) + PRINT *, 'tetamin = ', tetamin + idissip = max(1, int(tetamin / (2 * dtvr * iperiod))) * iperiod + PRINT *, 'idissip = ', idissip + dtdiss = idissip * dtvr + PRINT *, 'dtdiss = ', dtdiss, "s" + + call new_unit(unit) + open(unit, file="inidissip.csv", status="replace", action="write") + + ! Title line: + write(unit, fmt=*) '"presnivs (hPa)" "dtdiss * tetaudiv" ' & + // '"dtdiss * tetaurot" "dtdiss * tetah"' + + do l = 1, llm + write(unit, fmt=*) presnivs(l) / 100., dtdiss * tetaudiv(l), & + dtdiss * tetaurot(l), dtdiss * tetah(l) + end do + close(unit) + print *, 'Created file "inidissip.csv".' + + END SUBROUTINE inidissip + +end module inidissip_m