--- trunk/libf/dyn3d/Dissipation/inidissip.f90 2012/01/30 12:54:02 57 +++ trunk/libf/dyn3d/Dissipation/inidissip.f90 2012/09/20 09:57:03 65 @@ -8,7 +8,7 @@ REAL dtdiss ! in s integer idissip ! période de la dissipation (en pas de temps) - real tetaudiv(llm), tetaurot(llm), tetah(llm) ! in s + real tetaudiv(llm), tetaurot(llm), tetah(llm) ! in s-1 real cdivu, crot, cdivh contains @@ -20,16 +20,18 @@ ! Initialisation de la dissipation horizontale. Calcul des valeurs ! propres des opérateurs par méthode itérative. - USE comconst, ONLY : dtvr - use comdissnew, only: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, & - tetagrot, tetatemp - USE comvert, ONLY : preff, presnivs - USE conf_gcm_m, ONLY : iperiod - USE dimens_m, ONLY : iim, jjm, llm + USE comconst, ONLY: dtvr + use comdissnew, only: nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, & + tetatemp + USE comvert, 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 paramet_m, ONLY : jjp1 + use nxgraro2_m, only: nxgraro2 + USE paramet_m, ONLY: jjp1 ! Variables local to the procedure: REAL zvert(llm), max_zvert ! no dimension @@ -49,78 +51,51 @@ deltap = 1. call random_number(zh) zh = zh - 0.5 - CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE., 1) + CALL filtreg(zh, jjp1, 1, 2, 1, .TRUE.) DO l = 1, 50 - IF (lstardis) THEN - CALL divgrad2(1, zh, deltap, niterh, divgra, -1.) - ELSE - CALL divgrad(1, zh, niterh, divgra, -1.) - END IF - + CALL divgrad2(1, zh, deltap, niterh, divgra, -1.) zllm = abs(maxval(divgra)) zh = divgra / zllm END DO - IF (lstardis) THEN - cdivh = 1. / zllm - ELSE - cdivh = zllm**(- 1. / niterh) - END IF + 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., 1) + CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.) call random_number(zv) zv = zv - 0.5 - CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1) + CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.) DO l = 1, 50 - IF (lstardis) THEN - CALL gradiv2(zu, zv, nitergdiv, gx, gy, -1.) - ELSE - CALL gradiv(1, zu, zv, nitergdiv, gx, gy, -1.) - END IF - + CALL gradiv2(zu, zv, nitergdiv, gx, gy, -1.) zllm = max(abs(maxval(gx)), abs(maxval(gy))) zu = gx / zllm zv = gy / zllm end DO - IF (lstardis) THEN - cdivu = 1. / zllm - ELSE - cdivu = zllm**(- 1. / nitergdiv) - END IF + 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., 1) + CALL filtreg(zu, jjp1, 1, 2, 1, .TRUE.) call random_number(zv) zv = zv - 0.5 - CALL filtreg(zv, jjm, 1, 2, 1, .FALSE., 1) + CALL filtreg(zv, jjm, 1, 2, 1, .FALSE.) DO l = 1, 50 - IF (lstardis) THEN - CALL nxgraro2(1, zu, zv, nitergrot, gx, gy, -1.) - ELSE - CALL nxgrarot(1, zu, zv, nitergrot, gx, gy, -1.) - END IF - + CALL nxgraro2(zu, zv, nitergrot, gx, gy, -1.) zllm = max(abs(maxval(gx)), abs(maxval(gy))) zu = gx / zllm zv = gy / zllm end DO - IF (lstardis) THEN - crot = 1. / zllm - ELSE - crot = zllm**(- 1. / nitergrot) - END IF + crot = 1. / zllm PRINT *, 'crot = ', crot ! Variation verticale du coefficient de dissipation : @@ -131,15 +106,6 @@ tetaurot = zvert / tetagrot tetah = zvert / tetatemp - call new_unit(unit) - open(unit, file="inidissip.csv", status="replace", action="write") - write(unit, fmt=*) '"tetaudiv (s)" "tetaurot (s)" "tetah (s)"' ! title line - do l = 1, llm - write(unit, fmt=*) tetaudiv(l), tetaurot(l), tetah(l) - end do - close(unit) - print *, 'Created file "inidissip.csv".' - max_zvert = maxval(zvert) tetamin = min(1e6, tetagdiv / max_zvert, tetagrot / max_zvert, & tetatemp / max_zvert) @@ -149,6 +115,20 @@ 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