--- trunk/dyn3d/sortvarc.f 2014/03/05 14:57:53 82 +++ trunk/dyn3d/sortvarc.f 2014/09/04 10:05:52 104 @@ -2,64 +2,55 @@ IMPLICIT NONE + real, save:: ang, etot, ptot, ztot, stot, rmsdpdt, rmsv + contains - SUBROUTINE sortvarc(itau, ucov, teta, ps, masse, pk, phis, vorpot, phi, & - bern, dp, time_0) + SUBROUTINE sortvarc(ucov, teta, ps, masse, pk, phis, vorpot, phi, & + bern, dp, resetvarc) ! From dyn3d/sortvarc.F, version 1.1.1.1 2004/05/19 12:53:07 ! Author: P. Le Van ! Objet : sortie des variables de contrôle - USE comconst, ONLY: daysec, dtvr, g, omeg, rad - USE comgeom, ONLY: aire, cu, rlatu - USE conf_gcm_m, ONLY: day_step + USE comconst, ONLY: daysec, g, omeg, rad + USE comgeom, ONLY: aire_2d, cu_2d, rlatu USE dimens_m, ONLY: iim, jjm, llm - USE dynetat0_m, ONLY: day_ini - USE ener, ONLY: ang, ang0, etot, etot0, ptot, ptot0, rmsdpdt, rmsv, & - stot, stot0, ztot, ztot0 + USE ener, ONLY: ang0, etot0, ptot0, stot0, ztot0 use filtreg_m, only: filtreg use massbarxy_m, only: massbarxy - USE paramet_m, ONLY: iip1, iip2, ijp1llm, ip1jm, ip1jmp1, jjp1 + USE paramet_m, ONLY: iip1, ip1jm, jjp1 - INTEGER, INTENT(IN):: itau - REAL, INTENT(IN):: ucov(ip1jmp1, llm) - REAL, INTENT(IN):: teta(ip1jmp1, llm) - REAL, INTENT(IN):: ps(ip1jmp1) - REAL, INTENT(IN):: masse(ip1jmp1, llm) - REAL, INTENT(IN):: pk(ip1jmp1, llm) - REAL, INTENT(IN):: phis(ip1jmp1) + REAL, INTENT(IN):: ucov(iim + 1, jjm + 1, llm) + REAL, INTENT(IN):: teta(iim + 1, jjm + 1, llm) + REAL, INTENT(IN):: ps(iim + 1, jjm + 1) + REAL, INTENT(IN):: masse(iim + 1, jjm + 1, llm) + REAL, INTENT(IN):: pk(iim + 1, jjm + 1, llm) + REAL, INTENT(IN):: phis(iim + 1, jjm + 1) REAL, INTENT(IN):: vorpot(ip1jm, llm) - REAL, intent(in):: phi(ip1jmp1, llm) - real, intent(in):: bern(ip1jmp1, llm) - REAL, intent(in):: dp(ip1jmp1) - REAL, INTENT (IN):: time_0 + REAL, intent(in):: phi(iim + 1, jjm + 1, llm) + real, intent(in):: bern(iim + 1, jjm + 1, llm) + REAL, intent(in):: dp(iim + 1, jjm + 1) + logical, intent(in):: resetvarc ! Local: - REAL:: vor(ip1jm), bernf(ip1jmp1, llm), ztotl(llm) - REAL:: etotl(llm), stotl(llm), rmsvl(llm), angl(llm), ge(ip1jmp1) - REAL:: cosphi(ip1jm), omegcosp(ip1jm) - REAL dtvrs1j, rjour, heure, radsg, radomeg + REAL vor(ip1jm), bernf(iim + 1, jjm + 1, llm), ztotl(llm) + REAL etotl(llm), stotl(llm), rmsvl(llm), angl(llm), ge(iim + 1, jjm + 1) + REAL cosphi(2:jjm) + REAL radsg, radomeg REAL massebxy(ip1jm, llm) - INTEGER l, ij + INTEGER j, l, ij REAL ssum - real time !----------------------------------------------------------------------- PRINT *, "Call sequence information: sortvarc" - time = real(itau) / day_step + time_0 - dtvrs1j = dtvr/daysec - rjour = real(int(itau*dtvrs1j)) - heure = (itau*dtvrs1j-rjour)*24. - IF (abs(heure-24.)<=0.0001) heure = 0. - CALL massbarxy(masse, massebxy) ! Calcul de rmsdpdt ge = dp*dp - rmsdpdt = ssum(ip1jmp1, ge, 1) - ssum(jjp1, ge, iip1) + rmsdpdt = sum(ge) - sum(ge(1, :)) rmsdpdt = daysec*1.E-2*sqrt(rmsdpdt / (iim * jjp1)) bernf = bern CALL filtreg(bernf, jjp1, llm, -2, 2, .TRUE.) @@ -67,10 +58,7 @@ ! Calcul du moment angulaire radsg = rad/g radomeg = rad*omeg - DO ij = iip2, ip1jm - cosphi(ij) = cos(rlatu((ij-1)/iip1+1)) - omegcosp(ij) = radomeg*cosphi(ij) - END DO + cosphi = cos(rlatu(2:jjm)) ! Calcul de l'energie, de l'enstrophie, de l'entropie et de rmsv @@ -80,65 +68,51 @@ END DO ztotl(l) = (ssum(ip1jm, vor, 1)-ssum(jjm, vor, iip1)) - DO ij = 1, ip1jmp1 - ge(ij) = masse(ij, l) * (phis(ij) + teta(ij, l) * pk(ij, l) & - + bernf(ij, l)-phi(ij, l)) - END DO - etotl(l) = ssum(ip1jmp1, ge, 1) - ssum(jjp1, ge, iip1) - - DO ij = 1, ip1jmp1 - ge(ij) = masse(ij, l)*teta(ij, l) - END DO - stotl(l) = ssum(ip1jmp1, ge, 1) - ssum(jjp1, ge, iip1) - - DO ij = 1, ip1jmp1 - ge(ij) = masse(ij, l)*amax1(bernf(ij, l)-phi(ij, l), 0.) - END DO - rmsvl(l) = 2.*(ssum(ip1jmp1, ge, 1)-ssum(jjp1, ge, iip1)) - - DO ij = iip2, ip1jm - ge(ij) = (ucov(ij, l)/cu(ij)+omegcosp(ij))*masse(ij, l)*cosphi(ij) - END DO - angl(l) = radsg * (ssum(ip1jm-iip1, ge(iip2), 1) & - - ssum(jjm-1, ge(iip2), iip1)) + ge = masse(:, :, l) * (phis + teta(:, :, l) * pk(:, :, l) & + + bernf(:, :, l) - phi(:, :, l)) + etotl(l) = sum(ge) - sum(ge(1, :)) + + ge = masse(:, :, l)*teta(:, :, l) + stotl(l) = sum(ge) - sum(ge(1, :)) + + ge = masse(:, :, l) * max(bernf(:, :, l) - phi(:, :, l), 0.) + rmsvl(l) = 2.*(sum(ge)-sum(ge(1, :))) + + forall (j = 2:jjm) ge(:, j) = (ucov(:, j, l) / cu_2d(:, j) & + + radomeg * cosphi(j)) * masse(:, j, l) * cosphi(j) + angl(l) = radsg * (sum(ge(:, 2:jjm)) - sum(ge(1, 2:jjm))) END DO - DO ij = 1, ip1jmp1 - ge(ij) = ps(ij)*aire(ij) - END DO - ptot = ssum(ip1jmp1, ge, 1) - ssum(jjp1, ge, iip1) - etot = ssum(llm, etotl, 1) - ztot = ssum(llm, ztotl, 1) - stot = ssum(llm, stotl, 1) - rmsv = ssum(llm, rmsvl, 1) - ang = ssum(llm, angl, 1) - - IF (ptot0 == 0.) THEN - PRINT *, 'WARNING!!! On recalcule les valeurs initiales de :' - PRINT *, 'ptot, rmsdpdt, etot, ztot, stot, rmsv, ang' - PRINT *, ptot, rmsdpdt, etot, ztot, stot, rmsv, ang + ge = ps * aire_2d + ptot = sum(ge) - sum(ge(1, :)) + etot = sum(etotl) + ztot = sum(ztotl) + stot = sum(stotl) + rmsv = sum(rmsvl) + ang = sum(angl) + + IF (resetvarc .or. ptot0 == 0.) then + print *, 'sortvarc: recomputed initial values.' etot0 = etot ptot0 = ptot ztot0 = ztot stot0 = stot - ang0 = ang + ang0 = ang + PRINT *, 'ptot0 = ', ptot0 + PRINT *, 'etot0 = ', etot0 + PRINT *, 'ztot0 = ', ztot0 + PRINT *, 'stot0 = ', stot0 + PRINT *, 'ang0 = ', ang0 END IF - etot = etot/etot0 - rmsv = sqrt(rmsv/ptot) - ptot = ptot/ptot0 - ztot = ztot/ztot0 - stot = stot/stot0 - ang = ang/ang0 - - PRINT 3500, itau, int(day_ini + time), heure, time - PRINT 4000, ptot, rmsdpdt, etot, ztot, stot, rmsv, ang - -3500 FORMAT (4X, 'pas', I7, 5X, 'jour', i5, 1X, 'heure', F5.1, 4X, 'date', & - F10.5) -4000 FORMAT (10X, 'masse', 4X, 'rmsdpdt', 7X, 'energie', 2X, 'enstrophie', & - 2X, 'entropie', 3X, 'rmsv', 4X, 'mt.ang', /, 'GLOB ', F10.6, & - E13.6, 5F10.3/) + IF (.not. resetvarc) then + etot = etot/etot0 + rmsv = sqrt(rmsv/ptot) + ptot = ptot/ptot0 + ztot = ztot/ztot0 + stot = stot/stot0 + ang = ang/ang0 + end IF END SUBROUTINE sortvarc