Ignore:
Timestamp:
2021-07-27T16:46:06+02:00 (3 years ago)
Author:
agnes.ducharne
Message:

As done in r6136 for SPREAD_PREC_CONT, and in r6326 for interpolation of regional maps.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90

    r7257 r7258  
    154154  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)      :: time_precip 
    155155  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: timebnd_precip 
    156   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)      :: preciptime_slab             !! Variable needed to keep track of how much rainfall was already distributed 
     156  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: preciptime_slab             !! Variable needed to keep track of how much rainfall was already distributed 
    157157  ! 
    158158  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: swdown_slab, lwdown_slab 
    … …  
    925925    LOGICAL, SAVE :: first_call_spreadprec=.TRUE. 
    926926    REAL(r_std), SAVE :: time_to_spread=3600.0 
     927    LOGICAL, SAVE :: spreadprec_cont=.FALSE.  
     928    ! 
    927929    INTEGER(i_std) :: imin(1), i, tind(3) 
    928930    REAL(r_std) :: ft(3), dt, left, right 
    929931    INTEGER(i_std) :: offset, nb_spread 
    930932    LOGICAL, ALLOCATABLE, DIMENSION(:) :: mask 
     933    REAL(r_std), ALLOCATABLE, DIMENSION(:) :: tspread 
    931934    ! 
    932935    ! 
    … …  
    934937       ALLOCATE(mask(slab_size_max)) 
    935938       mask(:) = .FALSE. 
     939    ENDIF 
     940    IF ( .NOT. ALLOCATED(tspread) ) THEN  
     941       ALLOCATE(tspread(nbpoint_proc))  
     942       tspread(:) = time_to_spread 
    936943    ENDIF 
    937944    ! 
    … …  
    953960       ! 
    954961       IF ( nb_spread < 0 ) THEN 
    955           !Config Key   = SPRED_PREC_SEC 
     962          !Config Key   = SPREAD_PREC_SEC 
    956963          !Config Desc  = Spread the precipitation over an interval in seconds. 
    957964          !Config Def   = 3600 
    … …  
    967974          ! 
    968975          CALL getin_p('SPRED_PREC_SEC', time_to_spread) 
     976          IF ( time_to_spread == forcing_tstep_ave/2.0 ) THEN 
     977             ! If we still have the default value, check if the 
     978             ! keyword does not exist with the right spelling. 
     979             CALL getin_p('SPREAD_PREC_SEC', time_to_spread) 
     980          ENDIF 
    969981       ELSE 
    970982          time_to_spread = dt_sechiba_keep * nb_spread 
    971983       ENDIF 
     984       ! 
     985       ! Add the option to introduce a continuous distribution of precipitation 
     986       ! 
     987       !Config Key   = SPREAD_PREC_CONT 
     988       !Config Desc  = Take into account precipitation on the next forcing step for spreading it. 
     989       !Config Def   = FALSE 
     990       !Config Help  = This allows to extend the spreading of rainfall to the entire forcing 
     991       !Config         should it be raining on the following time step. This ensures that if it rains 
     992       !Config         in two consecutive forcing time step at least during the first period rainfall 
     993       !Config         will be continuous. This avoids the spiky nature of rainfall in periods of 
     994       !Config         prolonged rainfall. 
     995       !Config Units = - 
     996       ! 
     997       ! This is the default should 'SPREAD_PREC_CONT' not be present in the run.def 
     998       ! 
     999       spreadprec_cont = .FALSE. 
     1000       ! 
     1001       CALL getin_p('SPREAD_PREC_CONT', spreadprec_cont) 
    9721002       ! 
    9731003       ! Do some verifications on the information read from run.def 
    … …  
    10061036       ENDIF 
    10071037       ! 
    1008        ! 
    1009        ! 
    1010        ! Do we need to take some rain from the previous time step ? 
    1011        ! 
    1012        !! Time computation is not better than 1/1000 seconds 
    1013        IF ( time_int(1) < timebnd_central(tind(2),1) .AND. preciptime_slab(tind(1)) < (time_to_spread-0.001) ) THEN 
    1014           dt = ((timebnd_central(tind(2),1)-offset)-(time_int(1)-offset))*one_day 
    1015           ft(1) = MIN(time_to_spread - preciptime_slab(tind(1)), dt)/tlen 
    1016        ELSE 
    1017           ft(1) = zero 
    1018        ENDIF 
    1019        ! 
    1020        ! Is there still some rain to spread from the current forcing time step ? 
    1021        ! 
    1022        !! Time computation is not better than 1/1000 seconds 
    1023        IF (preciptime_slab(tind(2)) < (time_to_spread-0.001) ) THEN 
    1024           left = MAX(time_int(1), timebnd_central(tind(2),1)) 
    1025           right = MIN(time_int(2),timebnd_central(tind(2),2)) 
    1026           dt = ((right-offset)-(left-offset))*one_day 
    1027           ft(2) = MIN(time_to_spread - preciptime_slab(tind(2)), dt)/tlen 
    1028        ELSE 
    1029           ft(2) = zero 
    1030        ENDIF 
    1031        ! 
    1032        ! Do we need to take some rain from the next time step ? 
    1033        ! 
    1034        !! Time computation is not better than 1/1000 seconds 
    1035        IF ( time_int(2) > timebnd_central(tind(2),2) .AND. preciptime_slab(tind(3)) < (time_to_spread-0.001) ) THEN 
    1036           dt = ((time_int(2)-offset)-(timebnd_central(tind(2),2)-offset))*one_day 
    1037           ft(3) = MIN(time_to_spread - preciptime_slab(tind(3)), dt)/tlen 
    1038        ELSE 
    1039           ft(3) = zero 
    1040        ENDIF 
    1041        ! 
    1042        ! Do the actual calculation 
     1038       tspread(:) = time_to_spread 
    10431039       ! 
    10441040       DO i=1, nbpoint_proc 
     1041          ! 
     1042          IF ( spreadprec_cont ) THEN 
     1043             ! 
     1044             ! If the next time step has also rainfall, then time_to_spread becomes the length of the forcing time step. 
     1045             ! 
     1046             IF ( rainf_slab(i,tind(3)) > zero .OR. snowf_slab(i,tind(3)) > zero) THEN 
     1047                tspread(i) = forcing_tstep_ave 
     1048             ELSE 
     1049                tspread(i) = time_to_spread 
     1050             ENDIF 
     1051          ELSE 
     1052             ! Default behavious 
     1053             tspread(i) = time_to_spread 
     1054          ENDIF 
     1055          ! 
     1056          ! Do we need to take some rain from the previous time step ? 
     1057          ! 
     1058          !! Time computation is not better than 1/1000 seconds 
     1059          IF ( time_int(1) < timebnd_central(tind(2),1) .AND. preciptime_slab(i,tind(1)) < (tspread(i)-0.001) ) THEN 
     1060             dt = ((timebnd_central(tind(2),1)-offset)-(time_int(1)-offset))*one_day 
     1061             ft(1) = MIN(tspread(i) - preciptime_slab(i,tind(1)), dt)/tlen 
     1062          ELSE 
     1063             ft(1) = zero 
     1064          ENDIF 
     1065          ! 
     1066          ! Is there still some rain to spread from the current forcing time step ? 
     1067          ! 
     1068          !! Time computation is not better than 1/1000 seconds 
     1069          IF (preciptime_slab(i,tind(2)) < (tspread(i)-0.001) ) THEN 
     1070             left = MAX(time_int(1), timebnd_central(tind(2),1)) 
     1071             right = MIN(time_int(2),timebnd_central(tind(2),2)) 
     1072             dt = ((right-offset)-(left-offset))*one_day 
     1073             ft(2) = MIN(tspread(i) - preciptime_slab(i,tind(2)), dt)/tlen 
     1074          ELSE 
     1075             ft(2) = zero 
     1076          ENDIF 
     1077          ! 
     1078          ! Do we need to take some rain from the next time step ? 
     1079          ! 
     1080          !! Time computation is not better than 1/1000 seconds 
     1081          IF ( time_int(2) > timebnd_central(tind(2),2) .AND. preciptime_slab(i,tind(3)) < (tspread(i)-0.001) ) THEN 
     1082             dt = ((time_int(2)-offset)-(timebnd_central(tind(2),2)-offset))*one_day 
     1083             ft(3) = MIN(tspread(i) - preciptime_slab(i,tind(3)), dt)/tlen 
     1084          ELSE 
     1085             ft(3) = zero 
     1086          ENDIF 
     1087          ! 
     1088          ! Do the actual calculation 
     1089          ! 
    10451090          rainf(i) = (rainf_slab(i,tind(1)) * forcing_tstep_ave * ft(1) + & 
    1046                   &  rainf_slab(i,tind(2)) * forcing_tstep_ave * ft(2) + & 
    1047                   &  rainf_slab(i,tind(3)) * forcing_tstep_ave * ft(3))*tlen/time_to_spread 
    1048    
     1091               &  rainf_slab(i,tind(2)) * forcing_tstep_ave * ft(2) + & 
     1092               &  rainf_slab(i,tind(3)) * forcing_tstep_ave * ft(3))*tlen/tspread(i) 
     1093 
    10491094          snowf(i) = (snowf_slab(i,tind(1)) * forcing_tstep_ave * ft(1) + & 
    1050                   &  snowf_slab(i,tind(2)) * forcing_tstep_ave * ft(2) + & 
    1051                   &  snowf_slab(i,tind(3)) * forcing_tstep_ave * ft(3))*tlen/time_to_spread 
     1095               &  snowf_slab(i,tind(2)) * forcing_tstep_ave * ft(2) + & 
     1096               &  snowf_slab(i,tind(3)) * forcing_tstep_ave * ft(3))*tlen/tspread(i) 
     1097          ! 
     1098          ! Update the time over which we have already spread the rainf 
     1099          ! 
     1100          preciptime_slab(i,tind(1)) = preciptime_slab(i,tind(1)) + tlen * ft(1) 
     1101          preciptime_slab(i,tind(2)) = preciptime_slab(i,tind(2)) + tlen * ft(2) 
     1102          preciptime_slab(i,tind(3)) = preciptime_slab(i,tind(3)) + tlen * ft(3) 
     1103          ! 
    10521104       ENDDO 
    1053        ! 
    1054        ! Update the time over which we have already spread the rainf 
    1055        ! 
    1056        preciptime_slab(tind(1)) = preciptime_slab(tind(1)) + tlen * ft(1) 
    1057        preciptime_slab(tind(2)) = preciptime_slab(tind(2)) + tlen * ft(2) 
    1058        preciptime_slab(tind(3)) = preciptime_slab(tind(3)) + tlen * ft(3) 
    1059        ! 
    10601105    ELSE 
    10611106       WRITE(numout,*) "Time interval toward which we will interpolate : ", time_int 
    … …  
    12951340    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)    :: u_full, v_full 
    12961341    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)    :: ps_full, ztq_full, zuv_full 
     1342    REAL(r_std), ALLOCATABLE, DIMENSION(:)      :: preciptime_tmp 
    12971343    ! 
    12981344    ! 1.0 Verify that for the slabs the memory is allocated for the variable 
    … …  
    13151361    IF ( .NOT. ALLOCATED(time_precip) ) ALLOCATE(time_precip(slab_size)) 
    13161362    IF ( .NOT. ALLOCATED(timebnd_precip) ) ALLOCATE(timebnd_precip(slab_size,2)) 
    1317     IF ( .NOT. ALLOCATED(preciptime_slab) ) ALLOCATE(preciptime_slab(slab_size)) 
     1363    IF ( .NOT. ALLOCATED(preciptime_slab) ) ALLOCATE(preciptime_slab(nbpoint_proc,slab_size))  
     1364    IF ( .NOT. ALLOCATED(preciptime_tmp) ) ALLOCATE(preciptime_tmp(slab_size)) 
    13181365    ! 
    13191366    IF ( .NOT. ALLOCATED(swdown_slab) ) ALLOCATE(swdown_slab(nbpoint_proc,slab_size)) 
    … …  
    13641411            &                     tair_full, tairmax_full, tairmin_full, time_tair, timebnd_tair, & 
    13651412            &                     qair_full, time_qair, timebnd_qair, & 
    1366             &                     rainf_full, snowf_full, time_precip, timebnd_precip, preciptime_slab, & 
     1413            &                     rainf_full, snowf_full, time_precip, timebnd_precip, preciptime_tmp, & 
    13671414            &                     swdown_full, time_swdown, timebnd_swdown, & 
    13681415            &                     lwdown_full, time_lwdown, timebnd_lwdown, & 
    … …  
    14011448    CALL bcast(time_precip) 
    14021449    CALL bcast(timebnd_precip) 
    1403     CALL bcast(preciptime_slab) 
     1450    CALL bcast(preciptime_tmp) 
    14041451    CALL bcast(time_swdown) 
    14051452    CALL bcast(timebnd_swdown) 
    … …  
    14481495       zuv_slab(:,:) = zuv_full(:,:) 
    14491496    ENDIF 
     1497    ! 
     1498    ! 
     1499    ! 
     1500    DO is=1,nbpoint_proc 
     1501       preciptime_slab(is,:) = preciptime_tmp(:) 
     1502    ENDDO 
    14501503    ! 
    14511504    ! Clean-up to free the memory on the root processor. 
    … …  
    38613914        ! First extral the information from the date string 
    38623915        ! 
    3863         READ(ymd(1:INDEX(ymd,'-')-1),'(I4)') year0 
     3916        READ(ymd(1:INDEX(ymd,'-')-1),'(I)') year0 
    38643917        ymd=TRIM(ymd(INDEX(ymd,'-')+1:LEN_TRIM(ymd))) 
    3865         READ(ymd(1:INDEX(ymd,'-')-1),'(I2)') month0 
     3918        READ(ymd(1:INDEX(ymd,'-')-1),'(I)') month0 
    38663919        ymd=TRIM(ymd(INDEX(ymd,'-')+1:LEN_TRIM(ymd))) 
    3867         READ(ymd,'(I2)') day0 
     3920        READ(ymd,'(I)') day0 
    38683921        ! 
    38693922        ! Now extract from the time string 
    38703923        ! 
    3871         READ(hms(1:INDEX(hms,':')-1),'(I2)') hours0 
     3924        READ(hms(1:INDEX(hms,':')-1),'(I)') hours0 
    38723925        hms=TRIM(hms(INDEX(hms,':')+1:LEN_TRIM(hms))) 
    3873         READ(hms(1:INDEX(hms,':')-1),'(I2)') minutes0 
     3926        READ(hms(1:INDEX(hms,':')-1),'(I)') minutes0 
    38743927        hms=TRIM(hms(INDEX(hms,':')+1:LEN_TRIM(hms))) 
    3875         READ(hms,'(I2)') seci 
     3928        READ(hms,'(I)') seci 
    38763929        ! 
    38773930        sec0 = hours0*3600. + minutes0*60. + seci 
Note: See TracChangeset for help on using the changeset viewer.