Changeset 7262

Timestamp:

2021-07-27T19:02:48+02:00 (3 years ago)

Author:

agnes.ducharne

Message:

As in r7065 and r76158 for the new driver. Runs OK on jean-zay with the old driver, but not yet tested with the new driver. The corresponding configurations need to be added following r7081, r7082, and r7719.

Location:

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver

Files:

• forcing_tools.f90 (modified) (36 diffs)
• orchideedriver.f90 (modified) (10 diffs)

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90

@@ -204 +204 @@
     REAL(r_std) :: s_sec, e_sec, dateend, diff_sec, date_end
     INTEGER(i_std) :: i, ic
+    CHARACTER(LEN=20) :: str_cyclic_start(2), str_cyclic_end(2)
+    INTEGER(i_std) :: c_year_start, c_month_start, c_day_start, c_year_end, c_month_end, c_day_end
+    REAL(r_std) :: c_sec_start, c_sec_end
     !
     !Config Key  = START_DATE
@@ -212 +215 @@
     CALL getin('START_DATE',str_sdate)
     !
+    !Config Key  = CYCLIC_STARTDATE
+    !Config Desc = Date at which the cyclic year is started
+    !Config Def  = NONE
+    !Config Help = The format is the same as in the CF convention : 1999-09-13 12:0:0
+    str_cyclic_start = " "
+    CALL getin('CYCLIC_STARTDATE',str_cyclic_start)
+    
+    !
+    !Config Key  = CYCLIC_ENDDATE
+    !Config Desc = Date at which the cyclic year is ended
+    !Config Def  = NONE
+    !Config Help = The format is the same as in the CF convention : 1999-09-13 12:0:0
+    str_cyclic_end = " "
+    CALL getin('CYCLIC_ENDDATE',str_cyclic_end)
+    
+
+    !
+    ! the start date of simulation 
     IF ( (INDEX(str_sdate(1),"-") .NE. INDEX(str_sdate(1),"-", .TRUE.)) .AND. &
          &  (INDEX(str_sdate(2),":") .NE. INDEX(str_sdate(2),":", .TRUE.)) ) THEN
@@ -230 +251 @@
        CALL ipslerr(3, "forcing_integration_time", "START_DATE incorrectly specified in run.def", str_sdate(1), str_sdate(2))
     ENDIF
-    CALL ymds2ju (s_year, s_month, s_day, s_sec, date_start)
+    !---------------------------------
+    ! cyclic start date
+    IF ( (INDEX(str_cyclic_start(1),"-") .NE. INDEX(str_cyclic_start(1),"-", .TRUE.)) .AND. &
+         &  (INDEX(str_cyclic_start(2),":") .NE. INDEX(str_cyclic_start(2),":", .TRUE.)) ) THEN
+       DO i=1,2
+          tmpstr = str_cyclic_start(1)
+          ic = INDEX(tmpstr,"-")
+          tmpstr(ic:ic) = " "
+          str_cyclic_start(1) = tmpstr
+          tmpstr = str_cyclic_start(2)
+          ic = INDEX(tmpstr,":")
+          tmpstr(ic:ic) = " "
+          str_cyclic_start(2) = tmpstr
+       ENDDO
+       READ (str_cyclic_start(1),*) c_year_start, c_month_start, c_day_start
+       READ (str_cyclic_start(2),*) hours, minutes, seci
+       c_sec_start = hours*3600. + minutes*60. + seci
+    ELSE IF ( len_trim(str_cyclic_start(1)) .NE. 0 ) THEN
+       CALL ipslerr(3, "forcing_integration_time", "CYCLIC_STARTDATE incorrectly specified in run.def", str_cyclic_start(1), str_cyclic_start(2))
+    ENDIF
+    ! if s_year not the same as c_year, use cyear to compute date_start
+    IF ( ( s_year .NE. c_year_start) .AND. (len_trim(str_cyclic_start(1)) .NE. 0)) THEN
+       CALL ymds2ju (c_year_start, c_month_start, c_day_start, c_sec_start, date_start)
+    ELSE
+       CALL ymds2ju (s_year, s_month, s_day, s_sec, date_start)
+    ENDIF
     CALL forcing_printdate(date_start, "This is after reading the start date")
+
     !
     !Config Key  = END_DATE
@@ -258 +305 @@
        CALL ipslerr(3, "forcing_integration_time", "END_DATE incorrectly specified in run.def", str_edate(1), str_edate(2))
     ENDIF
-    CALL ymds2ju (e_year, e_month, e_day, e_sec, date_end)
-    !
-    CALL time_diff (s_year,s_month,s_day,s_sec,e_year,e_month,e_day,e_sec,diff_sec)
+
+    !---------------------------------
+    ! for cyclic end date
+    IF ( (INDEX(str_cyclic_end(1),"-") .NE. INDEX(str_cyclic_end(1),"-", .TRUE.)) .AND. &
+         &  (INDEX(str_cyclic_end(2),":") .NE. INDEX(str_cyclic_end(2),":", .TRUE.)) ) THEN
+       DO i=1,2
+          tmpstr = str_cyclic_end(1)
+          ic = INDEX(tmpstr,"-")
+          tmpstr(ic:ic) = " "
+          str_cyclic_end(1) = tmpstr
+          tmpstr = str_cyclic_end(2)
+          ic = INDEX(tmpstr,":")
+          tmpstr(ic:ic) = " "
+          str_cyclic_end(2) = tmpstr
+       ENDDO
+       READ (str_cyclic_end(1),*) c_year_end, c_month_end, c_day_end
+       READ (str_cyclic_end(2),*) hours, minutes, seci
+       c_sec_end = hours*3600. + minutes*60. + seci
+    ELSE IF ( len_trim(str_cyclic_end(1)) .NE. 0 ) THEN
+       CALL ipslerr(3, "forcing_integration_time", "CYCLIC_ENDDATE incorrectly specified in run.def", str_cyclic_end(1), str_cyclic_end(2))
+    ENDif
+
+    ! if e_year not the same as c_year_end, use cyear_end to compute date_end
+    IF (( e_year .NE. c_year_end)  .AND. (len_trim(str_cyclic_end(1)) .NE. 0) )THEN
+       CALL ymds2ju (c_year_end, c_month_end, c_day_end, c_sec_end, date_end)
+    ELSE
+       CALL ymds2ju (e_year, e_month, e_day, e_sec, date_end)
+    ENDIF
+    
+    !
+    IF (( s_year .NE. c_year_start) .AND. (len_trim(str_cyclic_start(1)) .NE. 0) )then
+       CALL time_diff (c_year_start,c_month_start,c_day_start,c_sec_start,c_year_end,c_month_end,c_day_end,c_sec_end,diff_sec)
+    ELSE
+       CALL time_diff (s_year,s_month,s_day,s_sec,e_year,e_month,e_day,e_sec,diff_sec)
+    ENDIF
+
     !
     !Config Key  = DT_SECHIBA
@@ -403 +483 @@
        ENDIF
        CALL forcing_getglogrid(nb_forcefile, forfilename, iim_tmp, jjm_tmp, nbpoint_tmp, .FALSE., landonly)
+       
        !
        IF ( PRESENT(wunit) ) THEN
@@ -429 +510 @@
           nbpoint_loc = nbpoint_in
        ENDIF
+       
        !
        ! Treat the time dimension now :
@@ -453 +535 @@
        ENDIF
        CALL forcing_integration_time(startdate, dt, nbdt)
-       !
+
        ! Test that the time interval requested by the user correspond to the time available in the 
        ! forcing file.
@@ -463 +545 @@
           CALL ipslerr (3,'forcing_open', 'Start time requested by the user is outside of the time interval',&
                & "covered by the forcing file.","Please verify the configuration in the run.def file.")
+          
        ENDIF
        !
@@ -496 +579 @@
        nbpoint_proc = nbpoint_glo
     ENDIF
+    
     !
     ! On the slave processes we need to allocate the memory for the data on root_prc to be bcast
@@ -748 +832 @@
     !
     mid_int = time_int(1) + (dt/2.0)/one_day
-    imin = MINLOC( ABS(time_central(1:slab_size) - mid_int), mask )
+    imin = MINLOC( ABS(time_central(1:slab_size) - mid_int), mask(1:slab_size) )
     !
     ! Verify that this is a possible date
@@ -757 +841 @@
        !
     ELSE
-       WRITE(*,*) "imin(1) = ", imin(1), (time_int_in(1) + (dt/2.0)/one_day)
        CALL forcing_printdate(time_int_in(1), "===> Start of target time interval.")
        CALL forcing_printdate(time_int_in(2), "===> End of target time interval.")
@@ -823 +906 @@
     mid_int = time_int(1) + (dt/2.0)/one_day
     ! Locate that time on the time axis of the forcing.
-    imin = MINLOC( ABS(time_central(1:slab_size) - mid_int), mask )
+    imin = MINLOC( ABS(time_central(1:slab_size) - mid_int), mask(1:slab_size) )
     !
     ! Determine which indices are to the left (slabind_a) and right (slabind_b) of the model time and will be used
     ! for the linear interpolation.
     !
+
     IF ( imin(1) > 1 .AND. imin(1) < slab_size ) THEN
        !
@@ -851 +935 @@
        IF ( mid_int < time_central(slabind_a) ) THEN
           IF ( time_int(2) < time_central(slabind_a) ) THEN
-             WRITE(*,*) "imin(1) = ", imin(1), (time_int_in(1) + (dt/2.0)/one_day)
              CALL forcing_printdate(time_int_in(1), "===> Start of target time interval.")
              CALL forcing_printdate(time_int_in(2), "===> End of target time interval.")
@@ -869 +952 @@
        IF ( mid_int > time_central(slabind_b) ) THEN
           IF ( time_int(1) > time_central(slabind_b) ) THEN
-             WRITE(*,*) "imin(1) = ", imin(1), (time_int_in(1) + (dt/2.0)/one_day)
              CALL forcing_printdate(time_int_in(1), "===> Start of target time interval.")
              CALL forcing_printdate(time_int_in(2), "===> End of target time interval.")
@@ -1027 +1109 @@
        ! In principle 3 time steps can contribute to the time step closest to the center of the forcing interval
        !
-       imin = MINLOC( ABS(time_central(1:slab_size)-(time_int(1)+time_int(2))/2.0), mask )
+       imin = MINLOC( ABS(time_central(1:slab_size)-(time_int(1)+time_int(2))/2.0), mask(1:slab_size) )
        tind(1) = MAX(imin(1)-1,1)
        tind(2) = imin(1)
@@ -1216 +1298 @@
     ! Locate the time step in the SLAB at hand
     !
-    imin = MINLOC( ABS(time_cent(1:slab_size)-(time_int(1)+time_int(2))/2.0), mask )
+    imin = MINLOC( ABS(time_cent(1:slab_size)-(time_int(1)+time_int(2))/2.0), mask(1:slab_size) )
     !
     ! Compute all the angels we will encounter for the current forcing interval
@@ -1241 +1323 @@
        julian_tmp = (time_int(1)+time_int(2))/2.0
        split_time = julian_tmp+split*tlen/one_day
-       tmin = MINLOC( ABS(time_cent(1:slab_size) - split_time), mask)
+       tmin = MINLOC( ABS(time_cent(1:slab_size) - split_time), mask(1:slab_size))
        DO WHILE (  tmin(1) .EQ. imin(1) .AND. split_time .LE. timebnd(slab_size,2) )
           split = split + 1
           split_time = julian_tmp+split*tlen/one_day
-          tmin = MINLOC( ABS(time_cent(1:slab_size) - split_time), mask)
+          tmin = MINLOC( ABS(time_cent(1:slab_size) - split_time), mask(1:slab_size))
        ENDDO
        !
@@ -1625 +1707 @@
        !
        first_call_readslab = .FALSE.
+       write(numout, *) "first_call_readslab in forcing_readslab_root"
        !
     ELSE       
@@ -1636 +1719 @@
        !
        current_offset = position_slab(2)-2
-       !
-    ENDIF
+       write(numout, *) "first_call_readslab in forcing_readslab_root 22"
+       !
+    ENDIF
+    
     !
     ! Check that the slab size is not too large
@@ -1653 +1738 @@
        inslabpos=1
        WRITE(*,*) ">> Reading from global position ", start_globtime, "up to ", end_globtime
+       write(numout,*) time_sourcefile 
        !
        DO if=MINVAL(time_sourcefile(start_globtime:end_globtime)),MAXVAL(time_sourcefile(start_globtime:end_globtime))
@@ -1693 +1779 @@
              CALL forcing_varforslab(if, "Tairmin", nctstart, nctcount, inslabpos, tairmin, cellmethod)
           ENDIF
+          
+          
           CALL forcing_varforslab(if, "Tair", nctstart, nctcount, inslabpos, tair, cellmethod)
           CALL forcing_attributetimeaxe(cellmethod, timeid_tair)
@@ -1796 +1884 @@
        !
        IF ( position_slab(2)-position_slab(1) .GT. slab_size ) THEN
-          WRITE(*,*) "Postition_slab =",  position_slab
-          WRITE(*,*) "Interval read : ", position_slab(2)-position_slab(1)
-          WRITE(*,*) "Time start and end : ", time(1,1), time(slab_size,1)
           DO it =1,nbtax
              WRITE(*,*) "Checking time_tmp on idex : ", it
@@ -1806 +1891 @@
              WRITE(*,*) "Interval read : ", imax(1)-imin(1)+1
           ENDDO
-          WRITE(*,*) "current_offset, slab_size =", current_offset, slab_size
           CALL ipslerr (3,'forcing_readslab_root',"The time slab read does not fit the number of variables read.",&
                &        "Could there be an error in the time axis ?"," ")
@@ -2411 +2495 @@
           !
        ENDDO
+       IF ( iim_glo == 0 .AND. jjm_glo == 0 ) THEN
+          CALL ipslerr (3,'forcing_getglogrid',"Did not recognize any dimensions in : ", filename(iff), " ")
+       ENDIF
     ENDDO
+    
     !
     ! 3.0 Read the spatial coordinate variables found in the first file.
@@ -2575 +2663 @@
        CALL forcing_contfrac2d(force_id(1), testvar_id, contfrac_id, testvar2d, contfrac2d, &
             & testvar_missing, contfrac_missing, nbland_glo)
+       
+
        !
        ! We have found a variable on which we can count the number of land points. We can build
@@ -2581 +2671 @@
        IF ( landonly ) THEN
           nbpoint_glo = nbland_glo
+          
+    
           IF ( .NOT. ALLOCATED(lindex_glo)) ALLOCATE(lindex_glo(nbpoint_glo))
           IF ( .NOT. ALLOCATED(contfrac_glo)) ALLOCATE(contfrac_glo(nbpoint_glo))
           IF ( .NOT. ALLOCATED(testvar)) ALLOCATE(testvar(nbpoint_glo))
+          
+ 
           IF ( contfrac_id > 0 ) THEN
              CALL forcing_buildindex(contfrac2d, contfrac_missing, landonly, lindex_glo, contfrac_glo)
              CALL forcing_reindex(iim_glo, jjm_glo, testvar2d, nbland_glo, testvar, lindex_glo)
+              
+
           ELSE
              CALL forcing_buildindex(testvar2d, testvar_missing, landonly, lindex_glo, testvar)
              contfrac_glo(:) = 1.0
+             
           ENDIF
        ELSE
@@ -2614 +2711 @@
     dump_mask = closefile 
     CALL forcing_checkindex(dump_mask, testvarname, testvar)
+    
     !
     !
@@ -2627 +2725 @@
     jjm_tmp = jjm_glo
     nbpoint_tmp = nbpoint_glo
+    
     !
     ! Clean up !
@@ -2668 +2767 @@
     !
     INTEGER(i_std) :: i,j,k
-    !
-    k=0
-    DO i=1,iim_glo
-       DO j=1,jjm_glo
-          IF ( landonly ) THEN
-             IF ( var2d(i,j) /= var_missing ) THEN
+
+    IF ( MAXVAL(var2d) >= var_missing ) THEN
+       ! Case when we have missing values in the vard2d
+       k=0
+       DO i=1,iim_glo
+          DO j=1,jjm_glo
+             IF ( landonly ) THEN
+                IF ( var2d(i,j) /= var_missing .AND. var2d(i,j) > 0.0 ) THEN
+                   k = k + 1
+                   lindex(k) = (j-1)*iim_glo+i 
+                   var(k) = var2d(i,j)
+                ENDIF
+             ELSE
+                ! When we take all point, no test is performed.
                 k = k + 1
                 lindex(k) = (j-1)*iim_glo+i 
                 var(k) = var2d(i,j)
              ENDIF
-          ELSE
-             ! When we take all point, no test is performed.
-             k = k + 1
-             lindex(k) = (j-1)*iim_glo+i 
-             var(k) = var2d(i,j)
-          ENDIF
+          ENDDO
        ENDDO
-    ENDDO
+    ELSE
+       ! We suppose that this is land fraction variable
+       k=0
+       DO i=1,iim_glo
+          DO j=1,jjm_glo
+             IF ( landonly ) THEN
+                IF ( var2d(i,j) > 0.0 ) THEN
+                   k = k + 1
+                   lindex(k) = (j-1)*iim_glo+i 
+                   var(k) = var2d(i,j)
+                ENDIF
+             ELSE
+                ! When we take all point, no test is performed.
+                k = k + 1
+                lindex(k) = (j-1)*iim_glo+i 
+                var(k) = var2d(i,j)
+             ENDIF
+          ENDDO
+       ENDDO
+    ENDIF
+    !
     !
   END SUBROUTINE forcing_buildindex
@@ -2721 +2843 @@
     ! First determine the contfrac variable
     !
+    
     IF ( contfrac_id > 0 ) THEN
        iret = NF90_INQUIRE_VARIABLE(ifile, contfrac_id, ndims=it)
@@ -2822 +2945 @@
                "We do not know how to handle this.", " ")
        ENDIF
-       !
-       ! Count the number of land points.
-       !
-       DO i=1,iim_glo
-          DO j=1,jjm_glo
-             IF ( contfrac(i,j) /= contfrac_missing ) THEN
-                nbland = nbland + 1
-             ENDIF
+       
+       IF ( MAXVAL(contfrac) >= contfrac_missing ) THEN
+          ! We have missing values in contfrac and we use it to count number of land points
+          DO i=1,iim_glo
+             DO j=1,jjm_glo
+                IF ( contfrac(i,j) /= contfrac_missing .AND. contfrac(i,j) > 0.0 ) THEN
+                   nbland = nbland + 1
+                ENDIF
+             ENDDO
           ENDDO
-       ENDDO
-       !
+          
+       ELSE
+          ! Then ocean is fully contfrc=0 !
+          DO i=1,iim_glo
+             DO j=1,jjm_glo
+                IF ( contfrac(i,j) > 0.0 ) THEN
+                   nbland = nbland + 1
+                ENDIF
+             ENDDO
+          ENDDO
+          
+       ENDIF
+
        ! If we did not find any land points on the map (i.e. iim_glo > 1 and jjm_glo > 1) then we
        ! look for fractions larger then 0.
@@ -2845 +2980 @@
           ENDDO
        ENDIF
-       !
+       
        ! Did we get a result ?
        !
@@ -2983 +3118 @@
        CALL forcing_reindex(nbpoint_glo, testvar, nbpoint_glo, testvar_reind, reindex_glo)
        !
+       
        CALL forcing_writetestvar("forcing_mask_glo.nc", iim_glo, jjm_glo, nbpoint_glo, &
             &                    lon_glo(:,1), lat_glo(1,:), lindex_glo, mask_glo, &
@@ -3248 +3384 @@
     !
     nbpoint_loc = SUM(mask_loc)
+    
+
     IF ( .NOT. zoom_forcing .AND. nbpoint_loc .NE. nbpoint_glo) THEN
        WRITE(*,*) "We have not zoomed into the forcing file still we get a number of"
@@ -3945 +4083 @@
         !
         DO it=1,nbtime_perfile(iff)
-           time_infiles(tstart+it) = date0_file(iff,tcnt) + time_read(it)*convtosec(iff)/one_day
+           !!time_infiles(tstart+it) = date0_file(iff,tcnt) + time_read(it)*convtosec(iff)/one_day
+           IF ( convtosec(iff) < one_day ) THEN
+              time_infiles(tstart+it) = date0_file(iff,tcnt) + time_read(it)*convtosec(iff)/one_day
+           ELSE
+              ! In the case of daily forcing the start time has to be 00UTC in Julian days.
+              time_infiles(tstart+it) = date0_file(iff,tcnt) + INT(time_read(it))
+           ENDIF
         ENDDO
         if ( check ) WRITE(*,*) "File ", iff, "goes from ",  time_infiles(tstart+1), " to ", &
@@ -4247 +4391 @@
         ! Zoom into the data and put it in the right place in the slab of data.
         !
+        
         CALL forcing_reindex(iim_glo, jjm_glo, timecount, tmp_slab2d, nbpoint_loc, slab_size, data, inslabpos, reindex2d_loc)
+        
      ENDIF
   ELSE

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/orchideedriver.f90

@@ -50 +50 @@
   USE ioipslctrl
   USE xios_orchidee
+  
   !
   !-
@@ -317 +318 @@
   IF ( is_root_prc) THEN
      CALL globgrd_getdomsz(gridfilename, iim_glo, jjm_glo, nbindex_g, model_guess, file_id, forfilename, zoom_lon, zoom_lat)
+     write(numout,*) "nbindex_g after calling globgrd_getdomsz in orchideedriver", nbindex_g
      nbseg = 4
   ENDIF
@@ -342 +344 @@
           &               kindex_g, contfrac_glo, calendar)
   ENDIF
+  
   !
   CALL bcast(lon_glo)
@@ -369 +372 @@
   ENDDO
   !
+  
   WRITE(*,*) "Rank", mpi_rank, " Before parallel region All land points : ",  nbindex_g
   WRITE(*,*) "Rank", mpi_rank, " from ", iim_glo, " point in Lon. and ", jjm_glo, "in Lat."
@@ -382 +386 @@
   !-
   CALL grid_set_glo(iim_glo, jjm_glo, nbindex_g)
+  
+
   CALL grid_allocate_glo(nbseg)
   ! Copy the list of indexes of land points into index_g used by ORCHIDEE and then broacast to all
@@ -477 +483 @@
   !-
   CALL forcing_integration_time(date0, dt, nbdt)
+  write(numout, *) "orchideedriver date0", date0, dt, nbdt
   !
   !-
@@ -558 +565 @@
   !
   CALL ioipslctrl_restini(itau, date0, dt, rest_id, rest_id_stom, itau_offset, date0_shifted)
-  WRITE(numout,*) "itau_offset : ", itau_offset, date0, date0_shifted
-  WRITE(numout,*) "itau_offset diff = ", date0_shifted, date0, date0_shifted-date0
   !
   ! To ensure that itau starts with 0 at date0 for the restart, we have to set an off-set to achieve this. 
@@ -574 +579 @@
        date0,    year_end, month_end, day_end, julian_diff,  &
        lon,      lat,      znt)
-
   CALL sechiba_xios_initialize 
 
@@ -599 +603 @@
      timestep_interval(1) = julian_start
      timestep_interval(2) = julian_end
-     julian = julian_end
+     julian = (julian_start + julian_end) /2.0  !julian_end
+     
      !
      ! Get the forcing data
@@ -605 +610 @@
      CALL forcing_getvalues(timestep_interval, dt, zlev_tq, zlev_uv, temp_air, qair, &
           &                 precip_rain, precip_snow, swdown, lwdown, sinang, u, v, pb)
+     
      !-
      !