Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90	(revision 7475)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90	(revision 7796)
@@ -51,4 +51,5 @@
   USE mod_orchidee_para
   USE forcingdaily_tools
+  USE grid
   !
   IMPLICIT NONE
@@ -90,5 +91,4 @@
   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,5 +404,5 @@
 !
   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,4 +424,5 @@
     INTEGER(i_std), OPTIONAL, INTENT(in)  :: wunit
     LOGICAL, OPTIONAL, INTENT(in)         :: landonly_arg
+    CHARACTER(LEN=*), INTENT(in)          :: model_guess
     !
     ! LOCAL
@@ -501,8 +502,8 @@
        !
        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,5 +1429,5 @@
     ! 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,5 +1436,8 @@
     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,5 +1461,7 @@
     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,4 +1488,13 @@
     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,8 +1521,8 @@
        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,5 +1561,4 @@
     CALL bcast(time_precip)
     CALL bcast(timebnd_precip)
-    CALL bcast(preciptime_tmp)
     CALL bcast(time_swdown)
     CALL bcast(timebnd_swdown)
@@ -1596,7 +1610,21 @@
     !
     !
-    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,8 +1658,8 @@
 !_ ==============================================================================================================================
 
-  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,4 +1673,5 @@
     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,5 +1686,4 @@
     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,6 +1724,4 @@
     CHARACTER(LEN=80) :: cellmethod
     !
-    LOGICAL, SAVE :: first_call_readslab=.TRUE.
-    !
     ALLOCATE(time_tmp(slab_size,nbtax))
     ALLOCATE(rau(nbpoint_loc,slab_size))
@@ -1713,5 +1739,4 @@
     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,18 +1747,10 @@
        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,5 +1947,4 @@
        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,5 +3269,5 @@
 !_ ==============================================================================================================================
 
-  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,5 +3278,5 @@
     !
     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,4 +3289,6 @@
     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,5 +3339,26 @@
     ! 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,4 +3385,10 @@
        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,21 +3428,45 @@
     ! 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,44 +3551,73 @@
     ! 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,5 +4109,4 @@
   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,7 +4359,4 @@
   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
