Changeset 1126 for LMDZ4/branches/LMDZ4-dev/libf/phylmd/geo2atm.F90

Timestamp:

Mar 17, 2009, 12:24:48 PM (15 years ago)

Author:

jghattas

Message:

• Correction de bug sur les poles pour le changement de repere avant et apres couplage avec l'ocean.
• Passage au F90

Olivier Marti + JG

File:

• LMDZ4/branches/LMDZ4-dev/libf/phylmd/geo2atm.F90 (modified) (2 diffs)

LMDZ4/branches/LMDZ4-dev/libf/phylmd/geo2atm.F90

@@ -5 +5 @@
   USE dimphy
   USE mod_phys_lmdz_para
-    
+
   IMPLICIT NONE
-  include 'dimensions.h'
+  INCLUDE 'dimensions.h'
+  INCLUDE 'YOMCST.h'
 
-! Change wind corrdinates from cartesian geocentric to local spherical
+! Change wind coordinates from cartesian geocentric to local spherical
 ! NB! Fonctionne probablement uniquement en MPI seul (sans OpenMP)
 !
@@ -17 +18 @@
   REAL, DIMENSION (im,jm), INTENT(OUT) :: pu, pv, pr
 
-  REAL, PARAMETER :: rpi = 3.141592653E0
-  REAL, PARAMETER :: rad = rpi / 180.0E0
+  REAL :: rad
+
+
+  rad = rpi / 180.0E0
   
-  REAL, DIMENSION (im,jm) :: zsinlon, zcoslon
-  REAL, DIMENSION (im,jm) :: zsinlat, zcoslat
+  pu(:,:) = &
+       - px(:,:) * SIN(rad * plon(:,:)) &
+       + py(:,:) * COS(rad * plon(:,:))
 
-  zsinlon = SIN (rad * plon)
-  zcoslon = COS (rad * plon)
-  zsinlat = SIN (rad * plat)
-  zcoslat = COS (rad * plat)
+  pv(:,:) = &
+       - px(:,:) * SIN(rad * plat(:,:)) * COS(rad * plon(:,:)) &
+       - py(:,:) * SIN(rad * plat(:,:)) * SIN(rad * plon(:,:)) &
+       + pz(:,:) * COS(rad * plat(:,:))  
 
-  pu = - px * zsinlon         + py * zcoslon
-  pv = - px * zsinlat*zcoslon - py * zsinlat*zsinlon + pz * zcoslat  
-  pr =   px * zcoslat*zcoslon + py * zcoslat*zsinlon + pz * zsinlat
+  pr(:,:) = &
+       + px(:,:) * COS(rad * plat(:,:)) * COS(rad * plon(:,:)) &
+       + py(:,:) * COS(rad * plat(:,:)) * SIN(rad * plon(:,:)) &
+       + pz(:,:) * SIN(rad * plat(:,:))
 
-! Value at North Pole
+  ! Value at North Pole
   IF (is_north_pole) THEN
-     pu(:,1) = - py(1,1)
-     pv(:,1) = - px(1,1) 
-     pr(:,1) = 0.
+     pu(:, 1) = pu(1, 1)
+     pv(:, 1) = pv(1, 1) 
+     pr(:, 1) = pr(1, 1)
   ENDIF
-
-! Value at South Pole     
+  
+  ! Value at South Pole     
   IF (is_south_pole) THEN
-     pu(:,jm) = py(1,jm)
-     pv(:,jm) = px(1,jm)
-     pr(:,jm) = 0.
+     pu(:,jm) = pu(1,jm)
+     pv(:,jm) = pv(1,jm)
+     pr(:,jm) = pr(1,jm)
   ENDIF
- 
+  
 END SUBROUTINE geo2atm