Changeset 8736

Timestamp:

2024-11-28T16:29:18+01:00 (3 hours ago)

Author:

jan.polcher

Message:

Backporting of Cassini projection and other improvements of the curvilinear grids to ORCHIDEE2.2. The code in the Trunk remains the reference.

Location:

branches/ORCHIDEE_2_2/ORCHIDEE

Files:

• src_driver/forcing_tools.f90 (modified) (1 diff)
• src_global/grid.f90 (modified) (10 diffs)

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90

@@ -3626 +3626 @@
              CALL grid_tolola(i_orig+0.5, j_orig-0.5, corners_loc(i,j,2,1), corners_loc(i,j,2,2)) 
              CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners_loc(i,j,3,1), corners_loc(i,j,3,2)) 
-             CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners_loc(i,j,4,1), corners_loc(i,j,4,2)) 
+             CALL grid_tolola(i_orig-0.5, j_orig+0.5, corners_loc(i,j,4,1), corners_loc(i,j,4,2)) 
              
              coslat = MAX( COS(lat_glo(i_orig,j_orig) * pi/180. ), mincos ) 

branches/ORCHIDEE_2_2/ORCHIDEE/src_global/grid.f90

@@ -1 @@
-
+!!
 !! This module define variables for the grid to gathered points.
 !!
-!! @call sechiba_main
-!! @Version : $Revision$, $Date$
-!! 
-!< $HeadURL$
-!< $Date$
-!< $Revision$
+!! $HeadURL$
+!! $Date$
+!! $Revision$
 !!
 !! @author Marie-Alice Foujols, Jan Polcher and Martial Mancip
@@ -31 +28 @@
   USE mod_orchidee_para_var
   USE mod_orchidee_transfert_para
+  
   USE haversine
   USE module_llxy
@@ -117 +115 @@
   !
   type (proj_info), SAVE, dimension(1:MAX_DOMAINS) :: proj_stack
+!$OMP THREADPRIVATE(proj_stack)
   !
   real(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)   :: dxwrf, dywrf
-  !
+!$OMP THREADPRIVATE(dxwrf, dywrf)
+  !
+  ! Full projected grid. It is needed to replace the information in the grid file which could
+  ! have been computed with anotehr Earth radius
+  !
+  REAL(r_std), ALLOCATABLE, DIMENSION (:,:), SAVE   :: lon_proj, lat_proj
+!$OMP THREADPRIVATE(lon_proj, lat_proj)
   !
   INTERFACE grid_tolola
@@ -742 +747 @@
     REAL(r_std)                :: dx, dy, dtx, dty, coslat
     REAL(r_std), ALLOCATABLE, DIMENSION (:)    :: LON, LAT
+    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)  :: LON_MAT, LAT_MAT
     REAL(r_std), ALLOCATABLE, DIMENSION (:,:)  :: mapfac_x, mapfac_y
+    REAL(r_std)                :: lon_test, lat_test, a, b
+    REAL(r_std)                :: maxerr=-10000000., minerr=10000000., d, xp, xg, yp, yg, zp, zg
+    LOGICAL                    :: checkgrid=.TRUE.
     !
     !
@@ -769 +778 @@
     iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DX", user_dxkm)
     iret = NF90_GET_ATT(fid, NF90_GLOBAL, "DY", user_dykm)
-    !
     user_dlat = user_dykm/(R_Earth * pi * 2.0 / 360.0)
     user_dlon = user_dxkm/(R_Earth * pi * 2.0 / 360.0)
@@ -775 +783 @@
     iret = NF90_GET_ATT(fid, NF90_GLOBAL, "POLE_LON", user_lon0)
     iret = NF90_GET_ATT(fid, NF90_GLOBAL, "POLE_LAT", user_lat0)
-
-    
+    !
+    ! Get one line or the full grid from file for test purposes
+    !
+    iret = NF90_inq_VARID(fid, "XLONG_M",lonid)
+    iret = NF90_inq_VARID(fid, "XLAT_M",latid)
+    iret = NF90_INQUIRE_VARIABLE(fid, lonid, dimids = dimIDs)
+    iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(1), len = numLons)
+    iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(2), len = numLats)
     !
     IF ( current_proj == PROJ_LATLON ) THEN
        !
-       iret = NF90_inq_VARID(fid, "XLONG_M",lonid)
-       iret = NF90_INQUIRE_VARIABLE(fid, lonid, dimids = dimIDs)
-       iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(1), len = numLons)
-       iret = NF90_INQUIRE_DIMENSION(fid, dimIDs(2), len = numLats)
        ALLOCATE(LON(numLons))
        iret = NF90_GET_VAR(fid, lonid, LON(:), start = (/ 1, 1, 1 /), count = (/ numLons, 1, 1 /))
-       
        iret = NF90_inq_VARID(fid, "XLAT_M",latid)
        ALLOCATE(LAT(numLats))
        iret = NF90_GET_VAR(fid, latid, LAT(:), start = (/ 1, 1, 1 /), count = (/ 1, numLats, 1 /))
-       
        user_dlon = (LON(numLons) - LON(1)) / (numLons - 1)
        user_dlat = (LAT(numLats) - LAT(1)) / (numLats - 1)
-       
-       DEALLOCATE(LON,LAT)
-
+    ELSE
+       !
+       ALLOCATE(LON_MAT(numLons,numLats))
+       iret = NF90_GET_VAR(fid, lonid, LON_MAT)
+       ALLOCATE(LAT_MAT(numLons,numLats))
+       iret = NF90_GET_VAR(fid, latid, LAT_MAT)
     ENDIF
     ! Unable to know from where to get the information
@@ -815 +826 @@
     ALLOCATE(dxwrf(iim,jjm))
     ALLOCATE(dywrf(iim,jjm))
+    ALLOCATE(lon_proj(iim,jjm))
+    ALLOCATE(lat_proj(iim,jjm)) 
     !
     varname = "MAPFAC_MX"
@@ -858 +871 @@
        
     else if (current_proj == PROJ_CASSINI) then
-       CALL rotate_coords(user_known_lat, user_known_lon, proj_known_lat, &
-            proj_known_lon, user_lat0, user_lon0, user_stand_lon, -1)
-       CALL map_set(current_proj, proj_stack(idom), &
-            loninc=user_dlon, &
-            latinc=user_dlat, &
-            lat1=user_known_lat, &
-            lon1=user_known_lon, &
-            lat0=user_lat0 ,&
-            lon0=user_lon0 ,&
-            knowni=user_known_x, &
-            knownj=user_known_y, &
-            stdlon=user_stand_lon, &
-            r_earth=R_Earth)
-
+                   CALL map_set(current_proj, proj_stack(idom),      &
+            latinc   = user_dlon,                        &
+            loninc   = user_dlat,                        &
+            lat1     = user_known_lat,                   &
+            lon1     = user_known_lon,                   &
+            lat0     = user_lat0,                        &
+            lon0     = user_lon0,                        &
+            knowni   = user_known_x,                     &
+            knownj   = user_known_y,                     &
+            stdlon   = user_stand_lon)
+       
     else if (current_proj == PROJ_LC) then
        call map_set(current_proj, proj_stack(idom), &
@@ -938 +948 @@
           !
           IF (proj_stack(idom)%code /= PROJ_LATLON ) THEN
-             dx = proj_stack(idom)%dx
-             ! Some projections in WRF do not store dy, in that case dy=dx.
-             IF ( proj_stack(idom)%dy > 0 ) THEN
-                dy = proj_stack(idom)%dy
+             IF ( proj_stack(idom)%dx > 0 ) THEN
+                dx = proj_stack(idom)%dx
+                ! Some projections in WRF do not store dy, in that case dy=dx.
+                IF ( proj_stack(idom)%dy > 0 ) THEN
+                   dy = proj_stack(idom)%dy
+                ELSE
+                   dy = proj_stack(idom)%dx
+                ENDIF
              ELSE
-                dy = proj_stack(idom)%dx
+                dx = user_dxkm
+                dy = user_dykm
              ENDIF
              dxwrf(i,j) = dx/mapfac_x(i,j)
@@ -963 +978 @@
     ENDDO
     !
+    ! Compute full projected grid
+    !
+    DO i=1,iim
+       DO j=1,jjm
+          CALL grid_tolola(i*1.0, j*1.0, lon_proj(i,j), lat_proj(i,j))
+       ENDDO
+    ENDDO
+    !
+    ! verify projection
+    !
+    IF ( checkgrid  .AND. current_proj .NE. PROJ_LATLON ) THEN
+       DO i=1,iim
+          DO j=1,jjm
+             d = haversine_distance(lon_proj(i,j), lat_proj(i,j), LON_MAT(i,j), LAT_MAT(i,j))
+             minerr=MIN(abs(d),minerr)
+             maxerr=MAX(abs(d),maxerr)
+          ENDDO
+       ENDDO
+       WRITE(*,*) "grid_initproj : Coordinate error range ", minerr, " < [m] < ", maxerr
+       IF ( maxerr > 1000 ) THEN
+          WRITE(*,*) "grid_initproj : The projected grid is different from the one found in grid.nc"
+          WRITE(*,*) "grid_initproj : This can be explained by a difference in Earth radius :"
+          WRITE(*,*) "grid_initproj : ORCHIDEE's Earth radius : ", R_Earth
+          WRITE(*,*) "grid_initproj : WRF's Earth radius (for instance) : 6370000."
+          !!          
+          !!  Assume that in this projection the central point of the grid
+          !!  is only affected by the radius of the Earth
+          !!
+          i=INT(iim/2.)
+          j=INT(jjm/2.)
+          !
+          CALL grid_tolola(i*1.0, j*1.0, lon_test, lat_test)
+          xp=COS(lat_test * pi/180.)*COS(lon_test * pi/180.)
+          xg=COS(LAT_MAT(i,j) * pi/180.)*COS(LON_MAT(i,j) * pi/180.)
+          yp=COS(lat_test * pi/180.)*SIN(lon_test * pi/180.)
+          yg=COS(LAT_MAT(i,j) * pi/180.)*SIN(LON_MAT(i,j) * pi/180.)
+          zp=SIN(lat_test * pi/180.)
+          zg=SIN(LAT_MAT(i,j) * pi/180.)
+          WRITE(*,*) "grid_initproj : Radius of grid file at mid-point : ", (R_Earth*xg/xp + R_Earth*yg/yp + R_Earth*zg/zp)/3.
+          !
+       ENDIF
+    ENDIF
+    !
+    !
+    !
+    IF ( ALLOCATED(LON) ) DEALLOCATE(LON)
+    IF ( ALLOCATED(LAT) ) DEALLOCATE(LAT)
+    IF ( ALLOCATED(LON_MAT) ) DEALLOCATE(LON_MAT)
+    IF ( ALLOCATED(LAT_MAT) ) DEALLOCATE(LAT_MAT)
+    !
   END SUBROUTINE grid_initproj
 !