wrfv2_fire/dyn_nmm/module_PRECIP_ADJUST.F

   1 !-----------------------------------------------------------------------
   2 !
   3 !NCEP_MESO:MODEL_LAYER: PHYSICS
   4 !
   5 !----------------------------------------------------------------------
   6 #include "nmm_loop_basemacros.h"
   7 #include "nmm_loop_macros.h"
   8 !-----------------------------------------------------------------------
   9 !
  10       MODULE MODULE_PRECIP_ADJUST
  11 !
  12 ! This module contains 3 subroutines:
  13 !     READPCP
  14 !     CHKSNOW
  15 !     ADJPPT
  16 !-----------------------------------------------------------------------
  17 !***
  18 !***  Specify the diagnostic point here: (i,j) and the processor number.
  19 !***  Remember that in WRF, local and global (i,j) are the same, so don't
  20 !***  use the "local(i,j)" output from glb2loc.f; use the GLOBAL (I,J)
  21 !***  and the PE_WRF.
  22 !***
  23 !
  24       INTEGER :: ITEST=346,JTEST=256,TESTPE=53
  25 !-----------------------------------------------------------------------
  26 !
  27       CONTAINS
  28 !
  29 !-----------------------------------------------------------------------
  30       SUBROUTINE READPCP(PPTDAT,DDATA,LSPA                              &
  31      &  ,IDS,IDE,JDS,JDE,KDS,KDE                                        &
  32      &  ,IMS,IME,JMS,JME,KMS,KME                                        &
  33      &  ,ITS,ITE,JTS,JTE,KTS,KTE)
  34 !
  35 !     ****************************************************************
  36 !     *                                                              *
  37 !     *   PRECIPITATION ASSIMILATION INITIALIZATION.                 *
  38 !     *    READ IN PRECIP ANALYSIS AND DATA MASK AND SET UP ALL      *
  39 !     *    APPROPRIATE VARIABLES.                                    *
  40 !     *                   MIKE BALDWIN, MARCH 1994                   *
  41 !     *                   Adapted to 2-D code, Ying Lin, Mar 1996    *
  42 !     *                   For WRF/NMM: Y.Lin Mar 2005                *
  43 !     *                                                              *
  44 !     ****************************************************************
  45 !-----------------------------------------------------------------------
  46 !
  47 ! READ THE BINARY VERSION OF THE PRECIP ANALYSIS.
  48 !
  49       IMPLICIT NONE
  50       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,                    &
  51      &                      IMS,IME,JMS,JME,KMS,KME,                    &
  52      &                      ITS,ITE,JTS,JTE,KTS,KTE
  53       REAL,DIMENSION(IDS:IDE,JDS:JDE) :: TEMPG
  54       REAL,DIMENSION(IMS:IME,JMS:JME) :: TEMPL
  55       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: DDATA, LSPA
  56       REAL,DIMENSION(IMS:IME,JMS:JME,3),INTENT(OUT) :: PPTDAT
  57       INTEGER :: I, J, IHR
  58       INTEGER :: MYPE
  59       CHARACTER*256 :: MESSAGE
  60 !
  61 ! Get the value of MYPE:
  62 !
  63       CALL WRF_GET_MYPROC(MYPE)
  64 !
  65       TEMPG=999.
  66 !
  67       DO IHR=1,3
  68         IF(MYPE==0)THEN
  69           READ(40+IHR) ((TEMPG(I,J),I=IDS,IDE-1),J=JDS,JDE-1)
  70           WRITE(MESSAGE,*) 'IHR=', IHR, ' FINISHED READING PCP TO TEMPG'
  71           CALL WRF_MESSAGE(MESSAGE)
  72           CLOSE(40+IHR)
  73 !
  74           DO J=JDS,JDE-1
  75             DO I=IDS,IDE-1
  76 ! In the binary version of the precip data, missing data are denoted as '999.'
  77 ! Convert the valid data from mm to m:
  78               IF (TEMPG(I,J).LT.900.) TEMPG(I,J)=TEMPG(I,J)*0.001
  79             ENDDO
  80           ENDDO
  81         ENDIF
  82 !
  83 ! Distribute to local temp array:
  84         CALL DSTRB(TEMPG,TEMPL,1,1,1,1,1                                &
  85      &,                IDS,IDE,JDS,JDE,KDS,KDE                          &
  86      &,                IMS,IME,JMS,JME,KMS,KME                          &
  87      &,                ITS,ITE,JTS,JTE,KTS,KTE)
  88 !
  89 ! Place into correct hour slot in PPTDAT:
  90         DO J=JMS,JME
  91           DO I=IMS,IME
  92             PPTDAT(I,J,IHR)=TEMPL(I,J)
  93           ENDDO
  94         ENDDO
  95 !
  96         IF(MYPE==TESTPE)THEN
  97           WRITE(MESSAGE,*) 'ADJPPT-READPCP, IHR',IHR, 'PPTDAT=',        &
  98      &      PPTDAT(ITEST,JTEST,IHR)
  99           CALL WRF_MESSAGE(MESSAGE)
 100         ENDIF
 101
 102       ENDDO
 103 !
 104 ! Give DDATA (hourly precipitation analysis partitioned into each physics
 105 ! timestep; partitioning done in ADJPPT) an initial value of 999, because
 106 ! TURBL/SURFCE is called before ADJPPT.  Also initialize LSPA to zero.
 107 !
 108       DDATA=999.
 109       LSPA=0.
 110 !
 111       RETURN
 112       END SUBROUTINE READPCP
 113 !
 114       SUBROUTINE CHKSNOW(NTSD,DT,NPHS,SR,PPTDAT                         &
 115      &  ,IDS,IDE,JDS,JDE,KDS,KDE                                        &
 116      &  ,IMS,IME,JMS,JME,KMS,KME                                        &
 117      &  ,ITS,ITE,JTS,JTE,KTS,KTE)
 118 !
 119 ! AT THE FIRST PHYSICS TIME STEP AFTER THE TOP OF EACH HOUR, CHECK THE SNOW
 120 ! ARRAY AGAINST THE SR (SNOW/TOTAL PRECIP RATIO).  IF SR .GE. 0.9, SET THIS
 121 ! POINT TO MISSING (SO WE WON'T DO SNOW ADJUSTMENT HERE).
 122 !
 123 !-----------------------------------------------------------------------
 124 !
 125       IMPLICIT NONE
 126 !
 127 !-----------------------------------------------------------------------
 128 !
 129       INTEGER,INTENT(IN) :: NTSD,NPHS
 130       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,                    &
 131      &                      IMS,IME,JMS,JME,KMS,KME,                    &
 132      &                      ITS,ITE,JTS,JTE,KTS,KTE
 133       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: SR
 134       REAL,DIMENSION(IMS:IME,JMS:JME,3),INTENT(INOUT) :: PPTDAT
 135       REAL,INTENT(IN) :: DT
 136       REAL :: TIMES
 137       INTEGER :: I, J, IHR
 138       INTEGER :: MYPE
 139       CHARACTER*256 :: MESSAGE
 140 !-----------------------------------------------------------------------
 141       TIMES=NTSD*DT
 142       IF (MOD(TIMES,3600.) < NPHS*DT) THEN
 143         IHR=INT(TIMES)/3600+1
 144         IF (IHR > 3) go to 10
 145         DO J=MYJS2,MYJE2
 146         DO I=MYIS1,MYIE1
 147           IF (SR(I,J) >= 0.9) PPTDAT(I,J,IHR) = 999.
 148         ENDDO
 149         ENDDO
 150 !
 151 ! Get the value of MYPE:
 152 !
 153         CALL WRF_GET_MYPROC(MYPE)
 154 !
 155         IF (MYPE==TESTPE) THEN
 156           WRITE(MESSAGE,1010) TIMES,SR(ITEST,JTEST)
 157  1010     FORMAT('ADJPPT-CHKSNOW: TIMES, SR=',F6.0,1X,F6.4)
 158           CALL WRF_MESSAGE(MESSAGE)
 159         ENDIF
 160       ENDIF
 161  10   CONTINUE
 162       RETURN
 163       END SUBROUTINE CHKSNOW
 164 !
 165       SUBROUTINE ADJPPT(NTSD,DT,NPHS,PREC,LSPA,PPTDAT,DDATA             &
 166      &  ,IDS,IDE,JDS,JDE,KDS,KDE                                        &
 167      &  ,IMS,IME,JMS,JME,KMS,KME                                        &
 168      &  ,ITS,ITE,JTS,JTE,KTS,KTE)
 169
 170 !***********************************************************************
 171 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
 172 !                .      .    .
 173 ! SUBPROGRAM:    ADJPPT    PRECIPITATION/CLOUD ADJUSTMENT
 174 !    PRGRMMR:    Y. LIN    ORG: W/NP22     DATE: 2005/03/30
 175 !
 176 ! ABSTRACT:
 177 !     ADJPPT  MAKES ADJUSTMENT TO MODEL'S TEMPERATURE, MOISTURE, HYDROMETEOR
 178 !     FIELDS TO BE MORE CONSISTENT WITH THE OBSERVED PRECIPITATION AND CLOUD
 179 !     TOP PRESSURE
 180 !
 181 !     FOR NOW, AS A FIRST STEP, JUST PARTITION THE INPUT HOURLY PRECIPITATION
 182 !     OBSERVATION INTO TIME STEPS, AND FEED IT INTO THE SOIL.
 183 ! PROGRAM HISTORY LOG:
 184 !
 185 !   2005/03/30  LIN      - BAREBONES PRECIPITATION PARTITION/FEEDING TO GROUND
 186 ! ATTRIBUTES:
 187 !   LANGUAGE: FORTRAN 90
 188 !   MACHINE : IBM
 189 !$$$
 190 !-----------------------------------------------------------------------
 191 !
 192       IMPLICIT NONE
 193 !
 194 !-----------------------------------------------------------------------
 195       INTEGER,INTENT(IN) :: NPHS, NTSD
 196       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,                    &
 197      &                      IMS,IME,JMS,JME,KMS,KME,                    &
 198      &                      ITS,ITE,JTS,JTE,KTS,KTE
 199       REAL,INTENT(IN) :: DT
 200       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: PREC
 201       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: DDATA, LSPA
 202       REAL,DIMENSION(IMS:IME,JMS:JME,3),INTENT(OUT) :: PPTDAT
 203 !-----------------------------------------------------------------------
 204 !***
 205 !***  LOCAL VARIABLES
 206 !***
 207 !-----------------------------------------------------------------------
 208       REAL :: DTPHS, FRACT, FRACT1, FRACT2, TIMES, TPHS1, TPHS2
 209       INTEGER :: I, J, IHR, IHR1, IHR2, NTSP
 210       INTEGER :: MYPE
 211       CHARACTER*256 :: MESSAGE
 212 !
 213 ! Get the value of MYPE:
 214 !
 215       CALL WRF_GET_MYPROC(MYPE)
 216 !
 217       TIMES=NTSD*DT
 218       IHR=INT(TIMES)/3600+1
 219 ! Size of physics time step:
 220       DTPHS=NPHS*DT
 221 !
 222 ! Compute the beginning and ending time of the current physics time step,
 223 ! TPHS1 and TPHS2:
 224 !
 225       NTSP=NTSD/NPHS+1
 226       TPHS1=(NTSP-1)*DTPHS
 227       TPHS2=NTSP*DTPHS
 228 !
 229       IHR1=INT(TPHS1)/3600+1
 230       IHR2=INT(TPHS2)/3600+1
 231 !
 232 ! Fraction of an hour that falls into IHR1 and IHR2.  Note that IHR1 and IHR2
 233 ! might be identical.
 234       IF (IHR1 > 3) THEN
 235         GO TO 200
 236       ELSEIF (IHR2 > 3) THEN
 237         IHR2=3
 238         FRACT1=(3600.- MOD(INT(TPHS1),3600))/3600.
 239         FRACT2=0.
 240       ELSEIF (IHR1 .EQ. IHR2) THEN
 241          FRACT1=0.5*DTPHS/3600.
 242          FRACT2=FRACT1
 243       ELSE
 244          FRACT1=(3600.- MOD(INT(TPHS1),3600))/3600.
 245          FRACT2=FLOAT(MOD(INT(TPHS2),3600))/3600.
 246       ENDIF
 247 !
 248       FRACT=FRACT1 + FRACT2
 249 !
 250       IF (MYPE==TESTPE) THEN
 251          WRITE(MESSAGE,1010) NTSD,NTSP,TIMES,IHR1,IHR2,TPHS1,TPHS2,      &
 252       &    FRACT1,FRACT2
 253  1010    FORMAT('ADJPPT: NTSD,NTSP,TIMES=',I4,1X,I4,1X,F6.0,' IHR1,IHR2=' &
 254       &   ,I1,1X,I1,' TPHS1,TPHS2=',F6.0,1X,F6.0,' FRACT1,FRACT2='        &
 255       &   ,2(1X,F6.4))
 256         CALL WRF_MESSAGE(MESSAGE)
 257       ENDIF
 258 !
 259 !-----------------------------------------------------------------------
 260 !   FRACT1/2 IS THE FRACTION OF IHR1/2'S PRECIP THAT WE WANT FOR
 261 !   THIS ADJUSTMENT (assuming that the physics time step spans over
 262 !   IHR1 and IHR2.  If not, then IHR1=IHR2).
 263 !-----------------------------------------------------------------------
 264 !   SET UP OBSERVED PRECIP FOR THIS TIMESTEP IN DDATA
 265 !-----------------------------------------------------------------------
 266       DO J=MYJS2,MYJE2
 267       DO I=MYIS1,MYIE1
 268 ! Note sometimes IHR1=IHR2.
 269         IF (PPTDAT(I,J,IHR1).GT.900..OR.PPTDAT(I,J,IHR2).GT.900.) THEN
 270           DDATA(I,J) = 999.
 271           LSPA(I,J) = LSPA(I,J) + PREC(I,J)
 272           GO TO 100
 273         ELSE
 274           IF (IHR2 .LE. 3) then
 275             DDATA(I,J) = PPTDAT(I,J,IHR1)*FRACT1                        &
 276      &                 + PPTDAT(I,J,IHR2)*FRACT2
 277           ELSE
 278             DDATA(I,J) = PPTDAT(I,J,IHR1)*FRACT1
 279           ENDIF
 280 !
 281            LSPA(I,J) = LSPA(I,J) + DDATA(I,J)
 282         ENDIF
 283         IF (I.EQ.ITEST .AND. J.EQ.JTEST .AND. MYPE.EQ.TESTPE) THEN
 284           WRITE(MESSAGE,1020) DDATA(I,J), PREC(I,J), LSPA(I,J)
 285  1020     FORMAT('ADJPPT: DDATA=',E12.6, ' PREC=',E12.6,' LSPA=',E12.6)
 286           CALL WRF_MESSAGE(MESSAGE)
 287         ENDIF
 288 !
 289  100    CONTINUE
 290       ENDDO
 291       ENDDO
 292 !
 293  200  CONTINUE
 294
 295       RETURN
 296       END SUBROUTINE ADJPPT
 297 END MODULE module_PRECIP_ADJUST