Changeset 7796 for branches/ORCHIDEE_2_2/ORCHIDEE/src_driver

Timestamp:

2022-11-08T16:24:46+01:00 (20 months ago)

Author:

xiaoni.wang

Message:

Updated the new driver in Tag2.2, with the corresponding changes in the commits 7714, 7770, 7771 and 7774 for trunk. Tested in debug and prod mode, with safran and cru forcings.

Location:

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver

Files:

• forcing_tools.f90 (modified) (28 diffs)
• globgrd.f90 (modified) (19 diffs)
• orchideedriver.f90 (modified) (10 diffs)

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90

@@ -51 +51 @@
   USE mod_orchidee_para
   USE forcingdaily_tools
+  USE grid
   !
   IMPLICIT NONE
@@ -90 +91 @@
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:)   :: time_bounds
   CHARACTER(LEN=20), SAVE, ALLOCATABLE, DIMENSION(:) :: time_axename, time_cellmethod
-  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)       :: preciptime
   INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:)    :: time_sourcefile
   INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:)  :: time_id
@@ -404 +404 @@
 !
   SUBROUTINE forcing_open(filenames_in, iim, jjm, lon, lat, nbpoint_in, drvzoom_lon, drvzoom_lat, &
-       &                  kindex, nbindex_perproc, wunit, landonly_arg)
+       &                  kindex, nbindex_perproc, wunit, model_guess, landonly_arg)
     !
     ! Opens the forcing file and reads some key information and stores them in the shared variables of the
@@ -424 +424 @@
     INTEGER(i_std), OPTIONAL, INTENT(in)  :: wunit
     LOGICAL, OPTIONAL, INTENT(in)         :: landonly_arg
+    CHARACTER(LEN=*), INTENT(in)          :: model_guess
     !
     ! LOCAL
@@ -501 +502 @@
        !
        IF ( PRESENT(wunit) ) THEN
-          WRITE(wunit,*) "Getting the zoomed grid", nbpoint_tmp
+          WRITE(wunit,*) "Getting the full grid", nbpoint_tmp
           CALL FLUSH(wunit)
        ENDIF
-       CALL forcing_zoomgrid(drvzoom_lon, drvzoom_lat, forfilename(1), .FALSE.)
+       CALL forcing_zoomgrid(drvzoom_lon, drvzoom_lat, forfilename(1), model_guess, .FALSE.)
        IF ( PRESENT(wunit) ) THEN
           WRITE(wunit,*) "Out of the zoomed grid operation"
@@ -1428 +1429 @@
     ! Local
     !
-    INTEGER(i_std)  :: is
+    INTEGER(i_std)  :: is, tlen
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)    :: tair_full, qair_full
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)    :: tairmin_full, tairmax_full
@@ -1435 +1436 @@
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)    :: u_full, v_full
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)    :: ps_full, ztq_full, zuv_full
-    REAL(r_std), ALLOCATABLE, DIMENSION(:)      :: preciptime_tmp
+    REAL(r_std), ALLOCATABLE, DIMENSION(:)      :: time_precip_old 
+    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)    :: savepreciptime 
+    INTEGER(i_std) :: s2in1(1) 
+    LOGICAL, SAVE :: first_call_readslab=.TRUE. 
     !
     ! 1.0 Verify that for the slabs the memory is allocated for the variable
@@ -1457 +1461 @@
     IF ( .NOT. ALLOCATED(timebnd_precip) ) ALLOCATE(timebnd_precip(slab_size,2))
     IF ( .NOT. ALLOCATED(preciptime_slab) ) ALLOCATE(preciptime_slab(nbpoint_proc,slab_size)) 
-    IF ( .NOT. ALLOCATED(preciptime_tmp) ) ALLOCATE(preciptime_tmp(slab_size))
+    IF ( .NOT. ALLOCATED(time_precip_old) ) ALLOCATE(time_precip_old(slab_size)) 
+    IF ( .NOT. ALLOCATED(savepreciptime) ) ALLOCATE(savepreciptime(nbpoint_proc,slab_size)) 
+    
     !
     IF ( .NOT. ALLOCATED(swdown_slab) ) ALLOCATE(swdown_slab(nbpoint_proc,slab_size))
@@ -1482 +1488 @@
     IF ( .NOT. ALLOCATED(zuv_slab) ) ALLOCATE(zuv_slab(nbpoint_proc,slab_size))
     !    
+    ! 
+    ! Save the old time axis for precipitation. 
+    ! 
+    IF ( first_call_readslab ) THEN 
+       time_precip_old(:) = zero 
+       preciptime_slab(:,:) = zero 
+    ELSE 
+       time_precip_old(:) = time_precip(:) 
+    ENDIF
     !
     IF ( is_root_prc) THEN
@@ -1506 +1521 @@
        ALLOCATE(zuv_full(nbpoint_loc,slab_size))
        !
-       CALL forcing_readslab_root(time_int, &
+       CALL forcing_readslab_root(time_int, first_call_readslab, &
             &                     tair_full, tairmax_full, tairmin_full, time_tair, timebnd_tair, &
             &                     qair_full, time_qair, timebnd_qair, &
-            &                     rainf_full, snowf_full, time_precip, timebnd_precip, preciptime_tmp, &
+            &                     rainf_full, snowf_full, time_precip, timebnd_precip, &
             &                     swdown_full, time_swdown, timebnd_swdown, &
             &                     lwdown_full, time_lwdown, timebnd_lwdown, &
@@ -1546 +1561 @@
     CALL bcast(time_precip)
     CALL bcast(timebnd_precip)
-    CALL bcast(preciptime_tmp)
     CALL bcast(time_swdown)
     CALL bcast(timebnd_swdown)
@@ -1596 +1610 @@
     !
     !
-    DO is=1,nbpoint_proc
-       preciptime_slab(is,:) = preciptime_tmp(:)
-    ENDDO
+    ! Transfer preciptime_slab from the old slab to the new one. 
+    ! s2in1: the index of the starting time from new slab in the previous slab 
+    s2in1 = MINLOC(ABS(time_precip_old-time_precip(1))) 
+    IF ( size(s2in1) .EQ. 0 ) THEN 
+       s2in1(1) = 1 
+    ENDIF
+    tlen=slab_size-s2in1(1)+1   
+    savepreciptime(:,:) = zero     
+    savepreciptime(:,1:tlen) = preciptime_slab(:,s2in1(1):slab_size)     
+    preciptime_slab(:,:) = zero 
+    preciptime_slab(:,1:tlen) = savepreciptime(:,1:tlen) 
+    
+    
+    IF (first_call_readslab) THEN 
+       first_call_readslab = .FALSE. 
+    ENDIF
+
     !
     ! Clean-up to free the memory on the root processor.
@@ -1630 +1658 @@
 !_ ==============================================================================================================================
 
-  SUBROUTINE forcing_readslab_root(time_int, &
+  SUBROUTINE forcing_readslab_root(time_int, first_call_readslab, &
             &                     tair, tairmax, tairmin, t_tair, tbnd_tair, &
             &                     qair, t_qair, tbnd_qair, &
-            &                     rainf, snowf, t_prec, tbnd_prec, prectime, &
+            &                     rainf, snowf, t_prec, tbnd_prec, &
             &                     swdown, t_swdown, tbnd_swdown, &
             &                     lwdown, t_lwdown, tbnd_lwdown, &
@@ -1645 +1673 @@
     REAL(r_std), INTENT(in)  :: time_int(2)                            !! The time interval over which the forcing is needed.
     !
+    LOGICAL,     INTENT(in)  :: first_call_readslab 
     REAL(r_std), INTENT(out) :: tair(:,:), tairmax(:,:), tairmin(:,:)
     REAL(r_std), INTENT(out) :: t_tair(:)
@@ -1657 +1686 @@
     REAL(r_std), INTENT(out) :: t_prec(:)
     REAL(r_std), INTENT(out) :: tbnd_prec(:,:)
-    REAL(r_std), INTENT(out) :: prectime(:)
     !
     REAL(r_std), INTENT(out) :: swdown(:,:)
@@ -1696 +1724 @@
     CHARACTER(LEN=80) :: cellmethod
     !
-    LOGICAL, SAVE :: first_call_readslab=.TRUE.
-    !
     ALLOCATE(time_tmp(slab_size,nbtax))
     ALLOCATE(rau(nbpoint_loc,slab_size))
@@ -1713 +1739 @@
     IF ( first_call_readslab ) THEN
        !
-       preciptime(:) = 0
        !
        ! If the first file is only there to provide the last time step of the previous year, we
@@ -1722 +1747 @@
        current_offset = MAX(imin(1)-2,1)
        !
-       first_call_readslab = .FALSE.
-       write(numout, *) "first_call_readslab in forcing_readslab_root"
        !
     ELSE       
        !
-       ! Put back the cummulated time of rainfall into the global array
-       !
-       preciptime(position_slab(1):position_slab(2)) = preciptime(position_slab(1):position_slab(2)) + &
-            &    prectime(1:slab_size)
-       !
        ! Compute new offset
        !
        current_offset = position_slab(2)-2
-       write(numout, *) "first_call_readslab in forcing_readslab_root 22"
        !
     ENDIF
@@ -1930 +1947 @@
        t_prec(1:slab_size) = time(position_slab(1):position_slab(2), timeid_precip)
        tbnd_prec(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_precip,:)
-       prectime(1:slab_size) = preciptime(position_slab(1):position_slab(2))
        !
        WRITE(*,*) "For SWdown we are using time axis '",TRIM(time_axename(timeid_swdown)),&
@@ -3253 +3269 @@
 !_ ==============================================================================================================================
 
-  SUBROUTINE forcing_zoomgrid (zoom_lon, zoom_lat, filename, dumpncdf)
+  SUBROUTINE forcing_zoomgrid (zoom_lon, zoom_lat, filename,model_guess, dumpncdf)
     !
     ! Get the area to be zoomed and the sizes of arrays we will need.
@@ -3262 +3278 @@
     !
     REAL(r_std), DIMENSION(2), INTENT(in) :: zoom_lon, zoom_lat
-    CHARACTER(LEN=*), INTENT(in) :: filename
+    CHARACTER(LEN=*), INTENT(in) :: filename,  model_guess 
     LOGICAL, INTENT(in) :: dumpncdf
     !
@@ -3273 +3289 @@
     REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_val, lat_val
     INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: zoom_index
+    INTEGER(i_std),DIMENSION(2) ::  tmp_lon_ind1, tmp_lat_ind1, tmp_lon_ind2, tmp_lat_ind2 
+    INTEGER(i_std) :: iindex_init, jindex_init, iindex_end, jindex_end, i_orig, j_orig 
     !
     INTEGER(i_std) :: iret, force_id, iv
@@ -3321 +3339 @@
     ! Determine the size in x and y of the zoom
     !
-    IF ( zoom_forcing ) THEN
+    lon_tmp(:) = 0.0 
+    lat_tmp(:) = 0.0 
+    IF ( zoom_forcing .AND. model_guess .EQ. 'WRF') THEN 
+       tmp_lon_ind1 = MINLOC(ABS(longlo_tmp(:,:)-zoom_lon(1))) 
+       tmp_lat_ind1 = MINLOC(ABS(lat_glo(:,:)-zoom_lat(1))) 
+       
+       tmp_lon_ind2 = MINLOC(ABS(longlo_tmp(:,:)-zoom_lon(2))) 
+       tmp_lat_ind2 = MINLOC(ABS(lat_glo(:,:)-zoom_lat(2))) 
+       
+       iindex_init = tmp_lon_ind1(1) 
+       jindex_init = tmp_lat_ind1(2) 
+                 
+       iindex_end= tmp_lon_ind2(1) 
+       jindex_end= tmp_lat_ind2(2) 
+       
+       lat_tmp(jindex_init:jindex_end) = 1.0 
+       lon_tmp(iindex_init:iindex_end) = 1.0 
+       
+       iim_loc = iindex_end - iindex_init +1 
+       jjm_loc = jindex_end - jindex_init +1 
+
+    ELSE IF ( zoom_forcing  .AND. model_guess .EQ. 'regular') THEN 
        !
        ! Working with the hypothesis it is a regular global grid and bring it back to the -180 to 180 interval
@@ -3346 +3385 @@
        jjm_loc = NINT(SUM(lat_tmp))
     ELSE
+       iindex_init = 1 
+       jindex_init = 1 
+       
+       iindex_end= size(longlo_tmp(:,1)) 
+       jindex_end= size(lat_glo(1,:)) 
+
        iim_loc = iim_glo
        jjm_loc = jjm_glo
@@ -3383 +3428 @@
     ! Now find the correspondance between the zoomed & re-ordered grid and the global one.
     !
-    DO i=1,iim_loc
-       DO j=1,jjm_loc
-          !
-          imin=MINLOC(ABS(longlo_tmp(:,1)-lon_val(i)))
-          jmin=MINLOC(ABS(lat_glo(1,:)-lat_val(j)))
-          !
-          lon_loc(i,j) = longlo_tmp(imin(1),jmin(1))
-          lat_loc(i,j) = lat_glo(imin(1),jmin(1))
-          mask_loc(i,j) = mask_glo(imin(1),jmin(1))
-          !
-          zoom_index(i,j,1) = imin(1)
-          zoom_index(i,j,2) = jmin(1)
-          !
+    IF (model_guess .EQ. 'WRF' .AND.  zoom_forcing) THEN 
+       ! if WRF grid and zoomed grid, added a few lines here 
+       mask_loc(:,:) = mask_glo(iindex_init:iindex_end, jindex_init:jindex_end) 
+       lat_loc(:,:) = lat_glo(iindex_init:iindex_end, jindex_init:jindex_end) 
+       lon_loc(:,:) = longlo_tmp(iindex_init:iindex_end, jindex_init:jindex_end) 
+       
+       
+       nbpoint_loc = 0 
+       DO i=1,iim_loc 
+          DO j=1,jjm_loc 
+             zoom_index(i,j,1) = iindex_init + i - 1 
+             zoom_index(i,j,2) = jindex_init + j - 1 
+             ! this is consistent to glogrd_getwrf code 
+             if ( mask_loc(i,j) .GT. 0.5 ) THEN 
+                nbpoint_loc = nbpoint_loc +1 
+             ENDIF
+          ENDDO
        ENDDO
-    ENDDO
-    !
-    nbpoint_loc = SUM(mask_loc)
+   
+   ELSE 
+      ! if regular grid, or wrf full grids: 
+      ! the orginal code, not modified 
+      DO i=1,iim_loc 
+         DO j=1,jjm_loc 
+            ! 
+            imin=MINLOC(ABS(longlo_tmp(:,1)-lon_val(i))) 
+            jmin=MINLOC(ABS(lat_glo(1,:)-lat_val(j))) 
+            ! 
+            lon_loc(i,j) = longlo_tmp(imin(1),jmin(1)) 
+            lat_loc(i,j) = lat_glo(imin(1),jmin(1)) 
+            mask_loc(i,j) = mask_glo(imin(1),jmin(1)) 
+            ! 
+            zoom_index(i,j,1) = imin(1) 
+            zoom_index(i,j,2) = jmin(1) 
+            ! 
+         ENDDO
+      ENDDO
+      ! 
+      nbpoint_loc = SUM(mask_loc) 
+      
+   ENDIF
     
 
@@ -3482 +3551 @@
     ! Treat first the longitudes
     !
-    DO j=1,jjm_loc
-       dx = zero
-       DO i=1,iim_loc-1
-          dx = dx+ABS(lon_loc(i,j)-lon_loc(i+1,j))
+    IF ( model_guess .EQ. "WRF" ) THEN 
+       ! if grid is wrf 
+       ! 
+       DO i=1,iim_loc 
+          DO j=1,jjm_loc 
+             i_orig = iindex_init + i - 1 
+             j_orig = jindex_init + j - 1 
+             
+             ! Corners 
+             CALL grid_tolola(i_orig+0.5, j_orig+0.5, corners_loc(i,j,1,1), corners_loc(i,j,1,2)) 
+             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)) 
+             
+             coslat = MAX( COS(lat_glo(i_orig,j_orig) * pi/180. ), mincos ) 
+             dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat 
+             dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth 
+             area_loc(i,j) = dx*dy 
+          ENDDO
        ENDDO
-       dx = dx/(iim_loc-1)
+
+    ELSE 
+       ! if regular grid: original code, not modified 
+       ! Treat first the longitudes 
+       ! 
+       DO j=1,jjm_loc 
+          dx = zero 
+          DO i=1,iim_loc-1 
+             dx = dx+ABS(lon_loc(i,j)-lon_loc(i+1,j)) 
+          ENDDO
+          dx = dx/(iim_loc-1) 
+          DO i=1,iim_loc 
+             corners_loc(i,j,1,1) = lon_loc(i,j)-dx/2.0 
+             corners_loc(i,j,2,1) = lon_loc(i,j)+dx/2.0 
+             corners_loc(i,j,3,1) = lon_loc(i,j)+dx/2.0 
+             corners_loc(i,j,4,1) = lon_loc(i,j)-dx/2.0 
+          ENDDO
+       ENDDO
+       ! 
+       ! Now treat the latitudes 
+       ! 
        DO i=1,iim_loc
-          corners_loc(i,j,1,1) = lon_loc(i,j)-dx/2.0
-          corners_loc(i,j,2,1) = lon_loc(i,j)+dx/2.0
-          corners_loc(i,j,3,1) = lon_loc(i,j)+dx/2.0
-          corners_loc(i,j,4,1) = lon_loc(i,j)-dx/2.0
+          dy = zero 
+          DO j=1,jjm_loc-1 
+             dy = dy + ABS(lat_loc(i,j)-lat_loc(i,j+1)) 
+          ENDDO
+          dy = dy/(jjm_loc-1) 
+          DO j=1,jjm_loc 
+             corners_loc(i,j,1,2) = lat_loc(i,j)+dy/2.0 
+             corners_loc(i,j,2,2) = lat_loc(i,j)+dy/2.0 
+             corners_loc(i,j,3,2) = lat_loc(i,j)-dy/2.0 
+             corners_loc(i,j,4,2) = lat_loc(i,j)-dy/2.0 
+          ENDDO
        ENDDO
-    ENDDO
-    !
-    ! Now treat the latitudes
-    !
-    DO i=1,iim_loc
-       dy = zero
-       DO j=1,jjm_loc-1
-          dy = dy + ABS(lat_loc(i,j)-lat_loc(i,j+1))
+    
+       ! 
+       ! Compute the areas of the grid (using the simplification that the grid is regular in lon/lat). 
+       ! 
+       DO i=1,iim_loc 
+          DO j=1,jjm_loc 
+             coslat = MAX( COS(lat_loc(i,j) * pi/180. ), mincos ) 
+             dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat 
+             dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth 
+             area_loc(i,j) = dx*dy 
+          ENDDO
        ENDDO
-       dy = dy/(jjm_loc-1)
-       DO j=1,jjm_loc
-          corners_loc(i,j,1,2) = lat_loc(i,j)+dy/2.0
-          corners_loc(i,j,2,2) = lat_loc(i,j)+dy/2.0
-          corners_loc(i,j,3,2) = lat_loc(i,j)-dy/2.0
-          corners_loc(i,j,4,2) = lat_loc(i,j)-dy/2.0
-       ENDDO
-    ENDDO
-    !
-    ! Compute the areas of the grid (using the simplification that the grid is regular in lon/lat).
-    !
-    DO i=1,iim_loc
-       DO j=1,jjm_loc
-          coslat = MAX( COS(lat_loc(i,j) * pi/180. ), mincos )
-          dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat
-          dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth
-          area_loc(i,j) = dx*dy
-       ENDDO
-    ENDDO
+       
+    ENDIF
+
     !
     ! If requested we read a variable, zoomin and dump the result into a test file.
@@ -4011 +4109 @@
   ALLOCATE(time(nb_forcing_steps, nbtax*nbtmethods), time_bounds(nb_forcing_steps,nbtax*nbtmethods,2))
   ALLOCATE(time_axename(nbtax*nbtmethods), time_cellmethod(nbtax*nbtmethods))
-  ALLOCATE(preciptime(nb_forcing_steps))
   ALLOCATE(time_sourcefile(nb_forcing_steps))
   ALLOCATE(time_id(nb_forcing_steps, nbtax))
@@ -4262 +4359 @@
   ENDDO
   !
-  ! Set to zero the variable for the cummulated time for rainfall
-  !
-  preciptime(:) = zero
   !
   ! Keep the very first date of the time axis for future use

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/globgrd.f90

@@ -80 +80 @@
     !
     INTEGER(i_std) :: iret, ndims, nvars, nb_atts, id_unlim, iv, lll
+    INTEGER(i_std) :: iindex_init, jindex_init, iindex_end, jindex_end
     INTEGER(i_std) :: iim_full, jjm_full, nbland_full
     CHARACTER(LEN=20) :: axname, varname
     CHARACTER(LEN=120) :: tmpfile
     REAL(r_std), DIMENSION(2) :: loczoom_lon, loczoom_lat
-    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask
+    INTEGER(i_std), DIMENSION(2) :: tmp_lon_ind1, tmp_lat_ind1, tmp_lon_ind2, tmp_lat_ind2 
+    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask,zoom_mask, lat, lon 
     !
     ! Set default values against which we can test
@@ -92 +94 @@
     !
     ! Verify the grid file name
+    ! if forcing_file exists and if grid file exists: tmpfile = grid file 
+    ! if forcing_file exits and grid file not exists: tmpfile = forcing_file 
+    ! if forcing_file not exists: tmp_file= grid file 
     ! 
-    IF ( INDEX(filename,"NONE") >= 1 ) THEN
+    IF ( PRESENT(forcingfile) ) THEN
        is_forcing_file=.TRUE.
-       IF ( PRESENT(forcingfile) ) THEN
+       IF ( INDEX(filename,"NONE") >= 1 ) THEN
           tmpfile=forcingfile(1)
        ELSE
-          CALL ipslerr (3,'globgrd_getdomsz',"Error No grid file provided :",tmpfile, " ")
+          tmpfile=filename
        ENDIF
     ELSE
@@ -105 +110 @@
     ENDIF
     !
-    ! Verify that the zomm is provided. Else choose the entire globe
+    ! Verify that the zoomed region is provided. Else choose the entire globe
     !
     IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN
@@ -144 +149 @@
           !! Coordinate variables used by WRF.
        CASE("west_east")
-          iim = lll
+          iim_full = lll
           model_guess = "WRF"
        CASE("south_north")
-          jjm = lll
+          jjm_full = lll
           model_guess = "WRF"
           !
           !! Variables used in WFDEI
        CASE("lon")
-          iim = lll
+          iim_full = lll
           model_guess = "regular"
        CASE("lat")
-          jjm = lll
+          jjm_full = lll
           model_guess = "regular"
        CASE("nbland")
-          nbland = lll
+          nbland_full = lll
           !
           !! Variables used by CRU-NCEP
        CASE("nav_lon")
-          iim = lll
+          iim_full = lll
           model_guess = "regular"
        CASE("nav_lat")
-          jjm = lll
+          jjm_full = lll
           model_guess = "regular"
        CASE("land")
-          nbland = lll
+          nbland_full = lll
+
+       
+          !! Variables used by CRU-JRA (v2.2.2) 
+       CASE("longitude") 
+          iim_full = lll 
+          model_guess = "regular" 
+       CASE("latitude") 
+          jjm_full = lll 
+          model_guess = "regular" 
        END SELECT
     ENDDO
     !
-    ! If we have a WRF file we need to count the number of land points.
+    ! If we have a WRF file we need to count the number of land points,  define iim jjm and mask 
     !
     IF (  model_guess == "WRF" ) THEN
 
-       ALLOCATE(mask(iim,jjm))
+       IF ( .NOT. ALLOCATED(mask) ) ALLOCATE(mask(iim_full,jjm_full)) 
 
        varname = "LANDMASK"
@@ -186 +200 @@
        ENDIF
 
-       nbland = COUNT(mask > 0.5)
-
-    ENDIF
-    !
-    !
-    ! If we are in the case of a forcing file a few functions from forcing_tools need to be called
+       nbland_full = COUNT(mask > 0.5) 
+       
+       
+       IF ( .NOT. ALLOCATED(lat)) ALLOCATE(lat(iim_full,jjm_full)) 
+       IF ( .NOT. ALLOCATED(lon)) ALLOCATE(lon(iim_full,jjm_full)) 
+       
+       varname = "XLONG_M" 
+       iret = NF90_INQ_VARID (fid, varname, iv) 
+       IF (iret /= NF90_NOERR) THEN 
+          CALL ipslerr (3,'globgrd_getdomsz',"Could not find variable ", varname," ") 
+       ELSE 
+          iret = NF90_GET_VAR (fid,iv,lon) 
+       ENDIF
+       
+       varname = "XLAT_M" 
+       iret = NF90_INQ_VARID (fid, varname, iv) 
+       IF (iret /= NF90_NOERR) THEN 
+          CALL ipslerr (3,'globgrd_getdomsz',"Could not find variable ", varname," ") 
+       ELSE 
+          iret = NF90_GET_VAR (fid,iv,lat) 
+       ENDIF
+       
+       !define nbland for full-region simulation in case of need 
+       nbland = nbland_full 
+       
+    ENDIF
+    ! 
+    ! 
+    ! If we are in the case of a forcing file and model_guess is regular grid, 
+    ! then a few functions from forcing_tools need to be called 
     ! so that the file is analysed with the tools of the forcing module.
     !
-    IF ( is_forcing_file ) THEN
+    ! predefine nbland, iim, jjm in the case is_forcing_file=false, or full grid simulation for wrf 
+    iim = iim_full 
+    jjm = jjm_full 
+    IF ( is_forcing_file .AND. model_guess .EQ. "regular") THEN 
        !
        ! Because we are re-using routines from the forcing module, we have to
@@ -204 +245 @@
        ENDIF
        !
-       ! This has to be a regular grid. A more clever clasification of files will be needed.
-       model_guess = "regular"
-
        ! Set last argument closefile=.FALSE. as the forcing file has been closed here above. 
        ! This will also induce that dump_mask=.FALSE. in forcing_getglogrid and the
@@ -213 +251 @@
        WRITE(*,*) forcingfile, "Forcing file with dimensions : ", iim_full, jjm_full, nbland_full
        !
-       CALL forcing_zoomgrid (loczoom_lon, loczoom_lat, forcingfile(1), .TRUE.)
+       CALL forcing_zoomgrid (loczoom_lon, loczoom_lat, forcingfile(1), model_guess, .TRUE.)  
        !
        CALL forcing_givegridsize (iim, jjm, nbland)
        !
+
+    ELSE IF ( is_forcing_file .AND. model_guess .EQ. "WRF") THEN 
+       ! 
+       ! if forcing_file exists and model_guess is wrf, and zoomed domain is asked, 
+       ! then we define the domain size according to the zoom 
+       ! if the full domain is asked, then we do not define the domain size again here. 
+       ! 
+       ! get the beginning and endding index in the case of zoomed region, for WRF grids (XW) 
+       ! Not to close the grid file here, to be used by globgrd_getgrid 
+       ! 
+       !IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN 
+       IF ( PRESENT(zoom_lon) .OR. PRESENT(zoom_lat) ) THEN 
+          
+          ! to get new iim, jjm and nbland for the zoomed region 
+          tmp_lon_ind1 = MINLOC(ABS(lon(:,:)-loczoom_lon(1))) 
+          tmp_lat_ind1 = MINLOC(ABS(lat(:,:)-loczoom_lat(1))) 
+          tmp_lon_ind2 = MINLOC(ABS(lon(:,:)-loczoom_lon(2))) 
+          tmp_lat_ind2 = MINLOC(ABS(lat(:,:)-loczoom_lat(2))) 
+          ! 
+          ! the zoomed region 
+          ! lon1, lat2------- lon2,lat2 
+          !   |                   | 
+          ! lon1, lat1------- lon2,lat1 
+          ! 
+          iindex_init = tmp_lon_ind1(1) 
+          jindex_init= tmp_lat_ind1(2) 
+          
+          iindex_end= tmp_lon_ind2(1) 
+          jindex_end= tmp_lat_ind2(2) 
+          
+          iim = iindex_end - iindex_init + 1 
+          jjm = jindex_end - jindex_init + 1 
+          
+          
+          IF ( .NOT. ALLOCATED(zoom_mask) ) ALLOCATE(zoom_mask(iim,jjm)) 
+          zoom_mask = mask(iindex_init:iindex_end, jindex_init:jindex_end) 
+          nbland = COUNT(zoom_mask > 0.5) 
+          IF ( ALLOCATED(zoom_mask) ) DEALLOCATE(zoom_mask) 
+          
+       ENDIF
+       
+       IF ( ALLOCATED(lat) ) DEALLOCATE(lat) 
+       IF ( ALLOCATED(lon) ) DEALLOCATE(lon) 
+       IF ( ALLOCATED(mask) ) DEALLOCATE(mask) 
+
     ENDIF
     !
@@ -250 +333 @@
   !---------------------------------------------------------------------
   SUBROUTINE globgrd_getgrid(fid, iim, jjm, nbland, model_guess, lon, lat, mask, area, corners, &
-       &                     lindex, contfrac, calendar)
+       &                     lindex, contfrac, calendar, zoom_lon, zoom_lat)
     !
     !
@@ -262 +345 @@
     INTEGER(i_std), INTENT(in)   :: iim, jjm, nbland
     CHARACTER(LEN=*), INTENT(in) :: model_guess
+    REAL(r_std), DIMENSION(2), INTENT(in), OPTIONAL :: zoom_lon, zoom_lat
     !
     ! OUTPUT
@@ -275 +359 @@
     CASE("WRF")
        CALL globgrd_getwrf(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
-            &               lindex, contfrac, calendar)
+            &               lindex, contfrac, calendar, zoom_lon, zoom_lat)
     CASE("regular")
        IF ( .NOT. is_forcing_file ) THEN
@@ -482 +566 @@
 !!
   SUBROUTINE globgrd_getwrf(fid, iim, jjm, nbland, lon, lat, mask, area, corners, &
-       &                     lindex, contfrac, calendar)
+       &                     lindex, contfrac, calendar, zoom_lon, zoom_lat)
     !
     USE defprec
@@ -491 +575 @@
     INTEGER(i_std), INTENT(in)   :: fid
     INTEGER(i_std), INTENT(in)   :: iim, jjm, nbland
+    REAL(r_std), DIMENSION(2), INTENT(in), OPTIONAL :: zoom_lon, zoom_lat
     !
     ! OUTPUT
@@ -503 +588 @@
     !
     INTEGER(i_std) :: i, ip, jp, k, iret, iv, nvars, varndim
-    INTEGER(i_std) :: imm1, imp1
+    INTEGER(i_std),DIMENSION(2) ::  tmp_lon_ind1, tmp_lat_ind1, tmp_lon_ind2, tmp_lat_ind2 
+    REAL(r_std), DIMENSION(2) :: loczoom_lon, loczoom_lat 
+    INTEGER(i_std) ::  ndims, nb_atts, id_unlim, lll, iim_full, jjm_full, iim_tmp, jjm_tmp 
+    INTEGER(i_std) :: iindex_init, jindex_init, iindex_end, jindex_end, i_orig, j_orig 
     CHARACTER(LEN=20) :: varname
-    REAL(r_std),DIMENSION(iim,jjm) :: mapfac_x, mapfac_y
+    
     INTEGER(i_std), DIMENSION(4) :: vardims
     INTEGER(i_std), DIMENSION(8) :: rose
     !
+    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_full, lon_full, lat_full, area_full 
+    REAL(r_std),ALLOCATABLE, DIMENSION(:,:) :: mapfac_x_full, mapfac_y_full 
+    CHARACTER(LEN=20) :: axname 
+    REAL(r_std) :: dx, dy, coslat 
+    ! 
+    REAL(r_std), PARAMETER :: mincos  = 0.0001 
+    REAL(r_std), PARAMETER :: pi = 3.141592653589793238 
+    REAL(r_std), PARAMETER :: R_Earth = 6378000. 
+    ! 
+    ! A lot of modifications are added to this subroutine (XW) 
     !
     ! Set some default values agains which we can check 
@@ -523 +621 @@
     calendar = 'gregorian'
     !
+    ! get dimension of full-grid data  
+    ! 
+    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars, & 
+         nAttributes=nb_atts, unlimitedDimId=id_unlim) 
+    ! 
+    DO iv=1,ndims 
+       ! 
+       iret = NF90_INQUIRE_DIMENSION (fid, iv, name=axname, len=lll) 
+       IF (iret /= NF90_NOERR) THEN 
+          CALL ipslerr (3,'globgrd_getdomsz',"Could not get size of dimension :"," "," ") 
+       ENDIF
+       ! 
+       ! This can be refined by testing the actual grid found in the file. 
+       ! 
+       SELECT CASE(axname) 
+          ! 
+          !! Coordinate variables used by WRF. 
+       CASE("west_east") 
+          iim_full = lll 
+          
+       CASE("south_north") 
+          jjm_full = lll 
+          
+       END SELECT
+    ENDDO
+
     !
     !  Init projection in grid.f90 so that it can be used later for projections.
     !
-    CALL grid_initproj(fid, iim, jjm)
+    CALL grid_initproj(fid, iim_full, jjm_full)
     !
     iret = NF90_INQUIRE (fid, nVariables=nvars)
@@ -533 +657 @@
        CALL ipslerr (3,'globgrd_getwrf',"Error inquiering variables from WRF grid file."," ", " ")
     ENDIF
+    IF ( .NOT. ALLOCATED(lon_full) ) ALLOCATE(lon_full(iim_full,jjm_full)) 
+    IF ( .NOT. ALLOCATED(lat_full) ) ALLOCATE(lat_full(iim_full,jjm_full)) 
+    IF ( .NOT. ALLOCATED(mask_full) ) ALLOCATE(mask_full(iim_full,jjm_full)) 
+    IF ( .NOT. ALLOCATED(mapfac_x_full) ) ALLOCATE(mapfac_x_full(iim_full,jjm_full)) 
+    IF ( .NOT. ALLOCATED(mapfac_y_full) ) ALLOCATE(mapfac_y_full(iim_full,jjm_full)) 
     !
     DO iv = 1,nvars
@@ -541 +670 @@
        !
        CASE("XLONG_M")
-          iret = NF90_GET_VAR(fid, iv, lon)
+          iret = NF90_GET_VAR(fid, iv, lon_full)
           IF (iret /= NF90_NOERR) THEN
              CALL ipslerr (3,'globgrd_getwrf',"Could not read the longitude from the WRF grid file.", " ", " ")
           ENDIF
        CASE("XLAT_M")
-          iret = NF90_GET_VAR(fid, iv, lat)
+          iret = NF90_GET_VAR(fid, iv, lat_full)
           IF (iret /= NF90_NOERR) THEN
              CALL ipslerr (3,'globgrd_getwrf',"Could not read the latitude from the WRF grid file.", " ", " ")
           ENDIF
        CASE("LANDMASK")
-          iret = NF90_GET_VAR(fid, iv, mask)
+          iret = NF90_GET_VAR(fid, iv, mask_full)
           IF (iret /= NF90_NOERR) THEN
              CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
           ENDIF
        CASE("MAPFAC_MX")
-          iret = NF90_GET_VAR(fid, iv, mapfac_x)
+          iret = NF90_GET_VAR(fid, iv, mapfac_x_full)
           IF (iret /= NF90_NOERR) THEN
              CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
           ENDIF
        CASE("MAPFAC_MY")
-          iret = NF90_GET_VAR(fid, iv, mapfac_y)
+          iret = NF90_GET_VAR(fid, iv, mapfac_y_full)
           IF (iret /= NF90_NOERR) THEN
              CALL ipslerr (3,'globgrd_getwrf',"Could not read the land mask from the WRF grid file.", " ", " ")
@@ -569 +698 @@
     ENDDO
     !
-    ! Compute corners and area on the iimxjjm grid
-    !
+    IF ( PRESENT(zoom_lon) .AND. PRESENT(zoom_lat) ) THEN  
+       !  
+       ! if zoomed region 
+       ! 
+       loczoom_lon = zoom_lon 
+       loczoom_lat = zoom_lat 
+       
+       
+       ! to get new iim, jjm and nbland for the zoomed region 
+       ! tmp_lon_ind1, tmp_lat_ind1 etc: dimension(2),  
+       ! with the first one for longitude, second one for latitude 
+       tmp_lon_ind1 = MINLOC(ABS(lon_full(:,:)-loczoom_lon(1))) 
+       tmp_lat_ind1 = MINLOC(ABS(lat_full(:,:)-loczoom_lat(1))) 
+       
+       tmp_lon_ind2 = MINLOC(ABS(lon_full(:,:)-loczoom_lon(2))) 
+       tmp_lat_ind2 = MINLOC(ABS(lat_full(:,:)-loczoom_lat(2))) 
+       ! 
+       ! get indices for the zoomed region 
+       ! lon1, lat2------- lon2,lat2 
+       !   |                   | 
+       ! lon1, lat1------- lon2,lat1 
+       ! 
+       iindex_init = tmp_lon_ind1(1) 
+       jindex_init = tmp_lat_ind1(2) 
+       
+       iindex_end= tmp_lon_ind2(1) 
+       jindex_end= tmp_lat_ind2(2) 
+       
+       ! allocate variable values for the zoomed region 
+       mask(:,:) = mask_full(iindex_init:iindex_end, jindex_init:jindex_end) 
+       lat(:,:) = lat_full(iindex_init:iindex_end, jindex_init:jindex_end) 
+       lon(:,:) = lon_full(iindex_init:iindex_end, jindex_init:jindex_end) 
+       
+       iim_tmp = iindex_end - iindex_init +1 
+       jjm_tmp = jindex_end - jindex_init +1 
+       
+       
+    ELSE 
+       ! 
+       ! if full grids 
+       ! 
+       iindex_init = 1 
+       jindex_init = 1  
+       
+       ! define variable values for full grids 
+       lat(:,:) = lat_full(:,:) 
+       lon(:,:) = lon_full(:,:) 
+       mask(:,:) = mask_full(:,:) 
+       
+    ENDIF                      
+    ! 
+    ! Compute corners on the iimxjjm full or zoomed grid 
     DO ip=1,iim
        DO jp=1,jjm
+          i_orig = iindex_init + ip - 1 
+          j_orig = jindex_init + jp - 1 
           ! Corners
-          CALL grid_tolola(ip+0.5, jp+0.5, corners(ip,jp,1,1), corners(ip,jp,1,2))
-          CALL grid_tolola(ip+0.5, jp-0.5, corners(ip,jp,2,1), corners(ip,jp,2,2))
-          CALL grid_tolola(ip-0.5, jp-0.5, corners(ip,jp,3,1), corners(ip,jp,3,2))
-          CALL grid_tolola(ip-0.5, jp-0.5, corners(ip,jp,4,1), corners(ip,jp,4,2))
+          CALL grid_tolola(i_orig+0.5, j_orig+0.5, corners(ip,jp,1,1), corners(ip,jp,1,2)) 
+          CALL grid_tolola(i_orig+0.5, j_orig-0.5, corners(ip,jp,2,1), corners(ip,jp,2,2)) 
+          CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners(ip,jp,3,1), corners(ip,jp,3,2)) 
+          CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners(ip,jp,4,1), corners(ip,jp,4,2)) 
           !
        ENDDO
     ENDDO
     !
-    ! Compute resolution and area on the gathered grid
+    ! Compute resolution and area on the gathered, full or zoomed grid 
     !
     k=0
@@ -588 +769 @@
     DO jp=1,jjm
        DO ip=1,iim
+          ! 
+          !get the right index of zoomed/full region in the original grids 
+          ! 
+          i_orig = iindex_init + ip - 1 
+          j_orig = jindex_init + jp - 1 
+          ! 
+          ! 
+          ! compute area 
+          coslat = MAX( COS(lat_full(i_orig,j_orig) * pi/180. ), mincos ) 
+          dx = ABS(corners(ip,jp,2,1) - corners(ip,jp,1,1)) * pi/180. * R_Earth * coslat 
+          dy = ABS(corners(ip,jp,1,2) - corners(ip,jp,3,2)) * pi/180. * R_Earth 
+          area(ip,jp) = dx*dy 
+          ! 
+          ! compute index and contfrac 
           IF ( mask(ip,jp) > 0.5 ) THEN
              !
+             ! index of the points in the local zoomed grid 
              k=k+1
              lindex(k) = (jp-1)*iim+ip 
@@ -602 +798 @@
        CALL ipslerr (3,'globgrd_getwrf',"Error closing the WRF grid file :", " ", " ")
     ENDIF
+    IF ( ALLOCATED(lon_full) ) DEALLOCATE(lon_full) 
+    IF ( ALLOCATED(lat_full) ) DEALLOCATE(lat_full) 
+    IF ( ALLOCATED(mask_full) ) DEALLOCATE(mask_full) 
+    IF ( ALLOCATED(mapfac_x_full) ) DEALLOCATE(mapfac_x_full) 
+    IF ( ALLOCATED(mapfac_y_full) ) DEALLOCATE(mapfac_y_full) 
+    
     !
   END SUBROUTINE globgrd_getwrf

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/orchideedriver.f90

@@ -26 +26 @@
 !================================================================================================================================
 ! 
-PROGRAM orchidedriver
+PROGRAM orchideedriver
   !---------------------------------------------------------------------
   !-
@@ -159 +159 @@
   ! Case when no ouput is desired.
   REAL(r_std), DIMENSION(2) :: testpt=(/9999.99,9999.99/)
-  INTEGER(i_std) :: ktest
+  INTEGER(i_std) :: ktest,alloc_stat 
   INTEGER :: printlev_loc  !! local write level
 
@@ -329 +329 @@
   !- Allocation of memory
   !- variables over the entire grid (thus in x,y)
-  ALLOCATE(lon_glo(iim_glo, jjm_glo))
-  ALLOCATE(lat_glo(iim_glo, jjm_glo))
-  ALLOCATE(mask_glo(iim_glo, jjm_glo))
-  ALLOCATE(area_glo(iim_glo, jjm_glo))
-  ALLOCATE(corners_glo(iim_glo, jjm_glo, nbseg, 2))
+  ALLOCATE(lon_glo(iim_glo, jjm_glo),stat=alloc_stat) 
+  ALLOCATE(lat_glo(iim_glo, jjm_glo),stat=alloc_stat) 
+  ALLOCATE(mask_glo(iim_glo, jjm_glo),stat=alloc_stat) 
+  ALLOCATE(area_glo(iim_glo, jjm_glo),stat=alloc_stat) 
+  ALLOCATE(corners_glo(iim_glo, jjm_glo, nbseg, 2),stat=alloc_stat) 
   !
   ! Gathered variables
-  ALLOCATE(kindex_g(nbindex_g))
-  ALLOCATE(contfrac_glo(nbindex_g))
-  !-
+  ALLOCATE(kindex_g(nbindex_g), stat=alloc_stat)
+  ALLOCATE(contfrac_glo(nbindex_g), stat=alloc_stat)
+  !--kindex_g: index of points in zoomed grids (if simulate zoomed region), or in full grids, both in case of regular grids 
   IF ( is_root_prc) THEN
      CALL globgrd_getgrid(file_id, iim_glo, jjm_glo, nbindex_g, model_guess, &
           &               lon_glo, lat_glo, mask_glo, area_glo, corners_glo,&
-          &               kindex_g, contfrac_glo, calendar)
+          &               kindex_g, contfrac_glo, calendar,  zoom_lon, zoom_lat)
   ENDIF
   
@@ -356 +356 @@
   CALL bcast(model_guess)
   !
-  ALLOCATE(lalo_glo(nbindex_g,2))
+  ALLOCATE(lalo_glo(nbindex_g,2), stat=alloc_stat)
   DO ik=1,nbindex_g
      !
@@ -386 +386 @@
   !-
   CALL grid_set_glo(iim_glo, jjm_glo, nbindex_g)
-  
+  CALL grid_allocate_glo(nbseg) 
 
-  CALL grid_allocate_glo(nbseg)
   ! Copy the list of indexes of land points into index_g used by ORCHIDEE and then broacast to all
   ! processors
@@ -397 +396 @@
   WRITE(numout,*) "Rank", mpi_rank, "Into Init_orchidee_data_para_driver with ", nbindex_g
   WRITE(numout,*) "Rank", mpi_rank, "Into ", index_g(1), index_g(nbindex_g)
+  CALL flush(numout) 
   !
   CALL Init_orchidee_data_para_driver(nbindex_g,index_g)
@@ -406 +406 @@
   ! Allocate grid on the local processor
   !
+  CALL flush(numout) 
+  ! 
+  ! Allocate grid varibles on the local processor by using nbp_loc:  
+  ! variables are: lalo, corners, area, neighbours, contfrac etc, with save attributes 
+  ! but these variables are not really used in the later parts of code. 
+  !  
   IF ( model_guess == "regular") THEN
      CALL grid_init (nbp_loc, nbseg, regular_lonlat, "ForcingGrid")
@@ -433 +439 @@
   !
   WRITE(numout,*) mpi_rank, "DIMENSIONS of grid on processor : iim, jjm, kjpindex = ", iim, jjm, kjpindex
+  CALL flush(numout) 
   !
   !  Allocate the local arrays we need :
@@ -498 +505 @@
   !
   CALL forcing_open(forfilename, iim_glo,  jjm_glo, lon_glo, lat_glo, nbindex_g, zoom_lon, zoom_lat, &
-       &            index_g, kjpindex, numout)
+       &            index_g, kjpindex, numout,  model_guess)
   !
   !-
@@ -972 +979 @@
   END SUBROUTINE orchideedriver_clear
   !
-END PROGRAM orchidedriver
+END PROGRAM orchideedriver