/[lmdze]/trunk/dyn3d/tourabs.f
ViewVC logotype

Diff of /trunk/dyn3d/tourabs.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/dyn3d/tourabs.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/dyn3d/tourabs.f90 revision 79 by guez, Fri Feb 28 17:52:47 2014 UTC
# Line 1  Line 1 
1        SUBROUTINE tourabs ( ntetaSTD,vcov, ucov, vorabs )  module tourabs_m
2    
3  c=======================================================================    IMPLICIT NONE
 c  
 c   Modif:  I. Musat (28/10/04)  
 c   -------  
 c   adaptation du code tourpot.F pour le calcul de la vorticite absolue  
 c   cf. P. Le Van  
 c  
 c   Objet:  
 c   ------  
 c  
 c    *******************************************************************  
 c    .............  calcul de la vorticite absolue     .................  
 c    *******************************************************************  
 c  
 c     ntetaSTD, vcov,ucov      sont des argum. d'entree pour le s-pg .  
 c             vorabs            est  un argum.de sortie pour le s-pg .  
 c  
 c=======================================================================  
   
       use dimens_m  
       use paramet_m  
       use comconst  
       use conf_gcm_m  
       use comgeom  
       use filtreg_m, only: filtreg  
       USE nr_util, ONLY : pi  
   
       IMPLICIT NONE  
 c  
       INTEGER ntetaSTD  
       REAL vcov( ip1jm,ntetaSTD ), ucov( ip1jmp1,ntetaSTD )  
       REAL vorabs( ip1jm,ntetaSTD )  
 c  
 c variables locales  
       INTEGER l, ij, i, j  
       REAL  rot( ip1jm,ntetaSTD )  
   
   
   
 c  ... vorabs = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) ..  
   
   
   
 c    ........  Calcul du rotationnel du vent V  puis filtrage  ........  
   
       DO 5 l = 1,ntetaSTD  
   
       DO 2 i = 1, iip1  
       DO 2 j = 1, jjm  
 c  
        ij=i+(j-1)*iip1  
        IF(ij.LE.ip1jm - 1) THEN  
 c  
         IF(cv(ij).EQ.0..OR.cv(ij+1).EQ.0..OR.  
      $     cu(ij).EQ.0..OR.cu(ij+iip1).EQ.0.) THEN  
          rot( ij,l ) = 0.  
          continue  
         ELSE  
          rot( ij,l ) = (vcov(ij+1,l)/cv(ij+1)-vcov(ij,l)/cv(ij))/  
      $                 (2.*pi*RAD*cos(rlatv(j)))*float(iim)  
      $                +(ucov(ij+iip1,l)/cu(ij+iip1)-ucov(ij,l)/cu(ij))/  
      $                 (pi*RAD)*(float(jjm)-1.)  
 c  
         ENDIF  
        ENDIF !(ij.LE.ip1jm - 1) THEN  
    2  CONTINUE  
   
 c    ....  correction pour  rot( iip1,j,l )  .....  
 c    ....     rot(iip1,j,l) = rot(1,j,l)    .....  
   
 CDIR$ IVDEP  
   
       DO 3 ij = iip1, ip1jm, iip1  
       rot( ij,l ) = rot( ij -iim, l )  
    3  CONTINUE  
   
    5  CONTINUE  
   
   
       CALL  filtreg( rot, jjm, ntetaSTD, 2, 1, .FALSE.)  
   
   
       DO 10 l = 1, ntetaSTD  
   
       DO 6 ij = 1, ip1jm - 1  
       vorabs( ij,l ) = ( rot(ij,l) + fext(ij)*unsairez(ij) )  
    6  CONTINUE  
   
 c    ..... correction pour  vorabs( iip1,j,l)  .....  
 c    ....   vorabs(iip1,j,l)= vorabs(1,j,l) ....  
 CDIR$ IVDEP  
       DO 8 ij = iip1, ip1jm, iip1  
       vorabs( ij,l ) = vorabs( ij -iim,l )  
    8  CONTINUE  
4    
5    10  CONTINUE  contains
6    
7        RETURN    SUBROUTINE tourabs(ntetaSTD, vcov, ucov, vorabs)
8        END  
9        ! Author: I. Musat, 28 October 2004, adaptation de tourpot
10        ! Objet : calcul de la vorticité absolue
11    
12        USE dimens_m, ONLY: iim, jjm
13        USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1
14        USE comconst, ONLY: rad
15        USE comgeom, ONLY: cu, cv, fext, rlatv, unsairez
16        USE filtreg_m, ONLY: filtreg
17        USE nr_util, ONLY: pi
18    
19        INTEGER, intent(in):: ntetaSTD
20        REAL, intent(in):: vcov(ip1jm, ntetaSTD), ucov(ip1jmp1, ntetaSTD)
21    
22        REAL, intent(out):: vorabs(ip1jm, ntetaSTD)
23        ! vorabs = Filtre(d(vcov)/dx - d(ucov)/dy) + fext
24    
25        ! Variables locales:
26        INTEGER l, ij, i, j
27        REAL rot(ip1jm, ntetaSTD)
28    
29        !--------------------------------------------------------------------
30    
31        ! Calcul du rotationnel du vent puis filtrage
32    
33        DO l = 1, ntetaSTD
34           DO i = 1, iip1
35              DO j = 1, jjm
36                 ij = i + (j - 1) * iip1
37                 IF (ij <= ip1jm - 1) THEN
38                    IF (cv(ij) == 0. .OR. cv(ij+1) == 0. .OR. cu(ij) == 0. &
39                         .OR. cu(ij+iip1) == 0.) THEN
40                       rot(ij, l) = 0.
41                    ELSE
42                       rot(ij, l) = (vcov(ij + 1, l) / cv(ij + 1) - vcov(ij, l) &
43                            / cv(ij)) / (2. * pi * RAD * cos(rlatv(j))) &
44                            * real(iim) + (ucov(ij + iip1, l) / cu(ij + iip1) &
45                            - ucov(ij, l) / cu(ij)) / (pi * RAD) * (real(jjm) - 1.)
46                    ENDIF
47                 ENDIF
48              end DO
49           end DO
50    
51           ! correction pour rot(iip1, j, l)
52           ! rot(iip1, j, l) = rot(1, j, l)
53           DO ij = iip1, ip1jm, iip1
54              rot(ij, l) = rot(ij - iim, l)
55           end DO
56        end DO
57    
58        CALL filtreg(rot, jjm, ntetaSTD, 2, 1, .FALSE.)
59    
60        DO l = 1, ntetaSTD
61           DO ij = 1, ip1jm - 1
62              vorabs(ij, l) = (rot(ij, l) + fext(ij) * unsairez(ij))
63           end DO
64    
65           ! correction pour vorabs(iip1, j, l)
66           ! vorabs(iip1, j, l)= vorabs(1, j, l)
67           DO ij = iip1, ip1jm, iip1
68              vorabs(ij, l) = vorabs(ij - iim, l)
69           end DO
70        end DO
71    
72      END SUBROUTINE tourabs
73    
74    end module tourabs_m
Legend:
Removed from v.76  
changed lines
  Added in v.79
  ViewVC Help
Powered by ViewVC 1.1.21