Forked from
Méso-NH / Méso-NH code
3259 commits behind the upstream repository.
-
WAUTELET Philippe authoredWAUTELET Philippe authored
write_ts1d.f90 13.76 KiB
!MNH_LIC Copyright 1995-2018 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
!! #############################
SUBROUTINE WRITE_TS1D
!! #############################
!! PURPOSE
!! -------
!! write quick and dirty 1D time series
!! METHOD
!! ------
!! write at each time step variables in PICASSO.EXE format
!! REFERENCE
!! ---------
!! none
!! AUTHOR
!! ------
!! PURPOSE
!! -------
!! write quick and dirty 1D time series
!! METHOD
!! ------
!! write at each time step variables in PICASSO.EXE format
!! REFERENCE
!! ---------
!! none
!! AUTHOR
!! ------
!! K. Suhre + C. Mari
!! MODIFICATIONS
!! -------------
!! Original 12/10/95
!! Add U, V, W, all water variables and all tracers to timeseries (KS)
!! 27/10/95 KS: add air density
!! 09/02/96 KS: change time to TDTSEG%TIME + ISCOUNT * XTSTEP
!! & grid identifier all set to 1 (mass levels)
!! 01/03/96 KS: remove air density
!! write variables every 900s to disk, regardless of XTSTEP
!! 02/03/96 KS: write time in hours rather than seconds
!! 24/04/96 KS: add PTIME to parameterlist
!! 29/07/96 KS: use XDUMMY1 variable from MODD_BLANK to define write-tstep
!! 31/07/96 KS: clip values that are smaller than 1E-40 to zero, since
!! otherwise strange behavior of PV-WAVE is possible
!! add printing of date and time and personalized comment
!! 07/08/96 KS: add printing of chemical species names and convert to ppt
!! print also XRHODREF so that conversiopn ppt to mol/cm3
!! could later be made
!! 20/02/97 KS: add cloud fraction XCLDFR
!! 05/03/97 KS: add writing of JNO2 and JO1D
!! 12/11/99 KS: add parallel open and close
!! 30/11/99 KS: add namelist parameters XCH_TS1D_TSTEP, CCH_TS1D_COMMENT,
!! and CCH_TS1D_FILENAME
!! 30/05/04 Tulet : update to be use with diag : now variables XCHEMLAT and
!! XCHEMLON could gives some profiles if specified in
!! namelist DIAG1.nam. Otherwise treatment is active only in
!! 1D case
!! 28/03/2018 P. Wautelet: replace TEMPORAL_DIST by DATETIME_DISTANCE
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!!!! EXTERNAL
!! --------
!! IMPLICIT ARGUMENTS
!! ------------------
USE MODE_DATETIME
USE MODE_FM, ONLY: IO_FILE_CLOSE_ll,IO_FILE_OPEN_ll
USE MODE_IO_MANAGE_STRUCT,ONLY: IO_FILE_ADD2LIST
USE MODE_IO_ll
USE MODE_GRIDPROJ
USE MODE_ll
!
USE MODD_LUNIT_n, ONLY: CINIFILE
USE MODD_NSV, ONLY: NSV,NSV_CHEMBEG,NSV_CHEMEND, &
NSV_AERBEG,NSV_AEREND
USE MODD_CH_AEROSOL, ONLY: CAERONAMES, LORILAM
USE MODD_DYN_n, ONLY: XTSTEP ! time-step of the model
USE MODD_DIM_n, ONLY: NKMAX ! # of points in Z of the physical grid
USE MODD_IO_ll, ONLY: TFILEDATA
USE MODD_PARAMETERS, ONLY: JPVEXT ! vertical external points number
USE MODD_GRID, ONLY: XLATORI,XLONORI
USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ
USE MODD_TIME_n, ONLY: TDTCUR ! Current Time and Date
USE MODD_TIME, ONLY: TDTEXP ! Begining Simulation Time and Date
USE MODD_CONF, ONLY: CPROGRAM,L1D ! Configuration of models
!!
USE MODD_FIELD_n, ONLY: XUT,XVT,XWT,&! wind speed at time t
XTHT, &! theta
XRT, &! water variables
XSVT, &! scalar variables
XTKET, &! Kinetic energy (rho e)
XCLDFR ! cloud fraction
USE MODD_CONF_n, ONLY: NRR ! Total number of water variables
USE MODD_REF_n, ONLY: XRHODREF ! dry air density of the reference state
USE MODD_CH_M9_n, ONLY: CNAMES ! names of chemical species
USE MODD_CH_MNHC_n, ONLY: LUSECHEM, &! flag that indicates if chemistry is used
XCH_TS1D_TSTEP, &! time between two call to write_ts1d
CCH_TS1D_COMMENT, &! comment for write_ts1d
CCH_TS1D_FILENAME ! filename for write_ts1d files
USE MODD_PARAM_n, ONLY: CRAD ! Kind of radiation parameterization
! 'NONE' if no parameterization
USE MODD_RADIATIONS_n, ONLY: XDTHRAD ! radiative heating/cooling rate
USE MODD_CH_JVALUES_n, ONLY: XJVALUES ! Jvalues and
USE MODD_CH_INIT_JVALUES, ONLY:JPJVMAX ! their number
USE MODD_PARAMETERS, ONLY: XUNDEF
USE MODD_DIAG_FLAG, ONLY: LCHEMDIAG, XCHEMLAT, XCHEMLON
IMPLICIT NONE
!! EXPLICIT ARGUMENTS
!! ------------------
! none
!! DECLARATION OF LOCAL VARIABLES
!! ------------------------------
REAL :: ZTIME
LOGICAL, SAVE :: GSFIRSTCALL = .TRUE.
INTEGER, SAVE :: ISIO1D ! IO-channel
CHARACTER(LEN=80), SAVE :: YSIO1DDEF ! name of def-file
CHARACTER(LEN=80), SAVE :: YSIO1DDAT ! name of dat-file
CHARACTER(LEN=40) :: YFORM = "(E15.8)" ! data output format
INTEGER, SAVE :: ISCOUNT = 0 ! timestep counter
INTEGER :: JK, JL, JN ! loop counter
INTEGER :: IIU,IJU ! domain size
INTEGER :: NBPROF ! number of profileINTEGER :: IUTIME ! hour in UT
INTEGER, SAVE :: ISSKIP ! # of timesteps to skip
INTEGER :: IINDEX, JINDEX ! index of each profile
REAL :: ZXINDEX, ZYINDEX ! distance from origin of each profile
REAL :: ZX, ZY ! poisition of each profile
REAL, DIMENSION(SIZE(XCHEMLAT)) :: ZLAT, ZLON
TYPE(TFILEDATA),POINTER,SAVE :: TZFILE
!
CHARACTER*8 :: YDATE ! for retrieval of date and time
CHARACTER*10 :: YTIME ! dito
CHARACTER*13 :: YLATLON ! dito
CHARACTER*13 :: YCLATLON ! dito
CHARACTER*4 :: YCYEAR ! current year
CHARACTER*2 :: YCMONTH ! current month
CHARACTER*2 :: YCDAY ! current day
CHARACTER*5 :: YCTIME ! current time
!!
!! EXECUTABLE STATEMENTS
!! ---------------------
!
TZFILE => NULL()
NBPROF = 0
CALL DATETIME_DISTANCE(TDTEXP,TDTCUR,ZTIME)
ZTIME = ZTIME + TDTEXP%TIME
CALL GET_DIM_EXT_ll ('B',IIU,IJU)
IF ((CPROGRAM =='DIAG ').AND.(LCHEMDIAG)) THEN
WRITE(YCYEAR,'(I4.4)') TDTCUR%TDATE%YEAR
WRITE(YCMONTH,'(I2.2)') TDTCUR%TDATE%MONTH
WRITE(YCDAY,'(I2.2)') TDTCUR%TDATE%DAY
IUTIME = INT(3600*MOD( 24.0+MOD((TDTCUR%TIME)/3600.,24.0),24.0 )) ! TU
WRITE(YCTIME,'(I5.5)') IUTIME
DO JN=1,SIZE(XCHEMLAT)
IF ((XCHEMLAT(JN) /= XUNDEF).AND.(XCHEMLON(JN) /= XUNDEF)) THEN
NBPROF = NBPROF + 1
ZLAT(NBPROF) = XCHEMLAT(JN)
ZLON(NBPROF) = XCHEMLON(JN)
END IF
END DO
END IF
IF (L1D) NBPROF = 1
DO JN=1,NBPROF
IF ((CPROGRAM =='DIAG ').AND.(LCHEMDIAG)) THEN
CALL SM_XYHAT(XLATORI,XLONORI, &
ZLAT(JN), ZLON(JN), ZX, ZY )
ZXINDEX=(ZX-XXHAT(1)) * IIU/(XXHAT(IIU)-XXHAT(1))
ZYINDEX=(ZY-XYHAT(1)) * IJU/(XYHAT(IJU)-XYHAT(1))
IINDEX=INT(ZXINDEX)
JINDEX=INT(ZYINDEX)
IF ((ABS(ZLAT(JN)) .GE. 10.).AND.(ABS(ZLON(JN)) .GE. 10.)) THEN
WRITE(YLATLON,'(F6.3,A1,F6.3)') ABS(ZLAT(JN)),'-',ABS(ZLON(JN))
WRITE(YCLATLON,'(A13)') TRIM(YLATLON)
END IF
IF ((ABS(ZLAT(JN)) .LT. 10.).AND.(ABS(ZLON(JN)) .GE. 10.)) THEN
WRITE(YLATLON,'(F5.3,A1,F6.3)') ABS(ZLAT(JN)),'-',ABS(ZLON(JN))
WRITE(YCLATLON,'(A12)') TRIM(YLATLON)
END IF
IF ((ABS(ZLAT(JN)) .LT. 10.).AND.(ABS(ZLON(JN)) .LT. 10.)) THEN
WRITE(YLATLON,'(F5.3,A1,F5.3)') ABS(ZLAT(JN)),'-',ABS(ZLON(JN))
WRITE(YCLATLON,'(A11)') TRIM(YLATLON)
END IF
IF ((ABS(ZLAT(JN)) .GE. 10.).AND.(ABS(ZLON(JN)) .LT. 10.)) THEN WRITE(YLATLON,'(F6.3,A1,F5.3)') ABS(ZLAT(JN)),'-',ABS(ZLON(JN))
WRITE(YCLATLON,'(A12)') TRIM(YLATLON)
END IF
CCH_TS1D_FILENAME = ADJUSTR(YCLATLON) // ":" // YCYEAR // YCMONTH // YCDAY // "-" // YCTIME
CCH_TS1D_COMMENT = "Fichier issu de " // CINIFILE
YSIO1DDEF = "def." // ADJUSTL(CCH_TS1D_FILENAME)
YSIO1DDAT = "dat." // ADJUSTL(CCH_TS1D_FILENAME)
END IF
IF (L1D) THEN
YSIO1DDEF = "def." // CCH_TS1D_FILENAME
YSIO1DDAT = "dat." // CCH_TS1D_FILENAME
IINDEX= 2
JINDEX= 2
END IF
IF ((IINDEX >= 1).AND.(IINDEX <= IIU).AND.&
(JINDEX >= 1).AND.(JINDEX <= IJU)) THEN
! write picasso def-file
IF (GSFIRSTCALL) THEN
CALL IO_FILE_ADD2LIST(TZFILE,YSIO1DDEF,'TXT','WRITE')
CALL IO_FILE_OPEN_ll(TZFILE,HPOSITION='REWIND',HSTATUS='NEW')
ISIO1D = TZFILE%NLU
! write comment
CALL DATE_AND_TIME(YDATE, YTIME)
WRITE(ISIO1D,*) YDATE, "@", YTIME, ": ", CCH_TS1D_COMMENT
! write # of points in Z
WRITE(ISIO1D,'(I5)') NKMAX
! write # of variables
WRITE(ISIO1D,'(I5)') 10 + NRR + NSV
! write grid identifier (all variables are assumed to be on mass points)
WRITE(ISIO1D,'(999A1)') ("1", JL = 1, 10 + NRR + NSV)
! write variable names
WRITE(ISIO1D,'(A)') "time (h)"
WRITE(ISIO1D,'(A)') "height (m)"
WRITE(ISIO1D,'(A)') "U (m/s)"
WRITE(ISIO1D,'(A)') "V (m/s)"
WRITE(ISIO1D,'(A)') "W (m/s)"
WRITE(ISIO1D,'(A)') "THT (K)"
WRITE(ISIO1D,'(A)') "TKE"
WRITE(ISIO1D,'(A)') "RHODREF (kg/m3)"
WRITE(ISIO1D,'(A)') "XDTHRAD radiat. heating/cooling rate"
WRITE(ISIO1D,'(A)') "XCLDFR (cloud fraction)"
WRITE(ISIO1D,'(A)') "JNO2 (1/s)"
WRITE(ISIO1D,'(A)') "JO1D (1/s)"
DO JL = 1, NRR
WRITE(ISIO1D,'(A, I1, A)') "Q", JL, " water variable (kg/kg)"
ENDDO
DO JL=1,NSV
IF (LUSECHEM .AND. JL >= NSV_CHEMBEG .AND. JL <= NSV_CHEMEND) THEN
WRITE(ISIO1D,'(A)') TRIM(CNAMES(JL-NSV_CHEMBEG+1)) // " (ppt)"
ELSE IF(LORILAM .AND. JL >= NSV_AERBEG .AND. JL <= NSV_AEREND) THEN
WRITE(ISIO1D,'(A)') TRIM(CAERONAMES(JL-NSV_AERBEG+1)) // " (ppt)"
ELSE
WRITE(ISIO1D,'(A, I3.3, A)') "CHEM", JL, " scalar variable (kg/kg)"
END IF
END DO
CALL IO_FILE_CLOSE_ll(TZFILE)
TZFILE => NULL()
! open picasso dat-file CALL IO_FILE_ADD2LIST(TZFILE,YSIO1DDAT,'TXT','WRITE')
CALL IO_FILE_OPEN_ll(TZFILE,HPOSITION='REWIND',HSTATUS='NEW')
ISIO1D = TZFILE%NLU
! calculate ISSKIP
IF (L1D) THEN
ISSKIP = MAX(1,IFIX(XCH_TS1D_TSTEP/XTSTEP))
END IF
IF ((CPROGRAM =='DIAG ').AND.(LCHEMDIAG)) ISSKIP = 1
PRINT *, "WRITE_TS1D: every ", ISSKIP, " time steps will be stored (!)"
END IF
END IF
IF (MOD(ISCOUNT,ISSKIP).EQ.0) THEN
! write a single data set at time t (in hours)
CALL WRITECLIP ( ZTIME / 3600. )
! print height variable
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XZZ(IINDEX,JINDEX,JK) )
ENDDO
! print U
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XUT(IINDEX,JINDEX,JK) )
ENDDO
! print V
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XVT(IINDEX,JINDEX,JK) )
ENDDO
! print W
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XWT(IINDEX,JINDEX,JK) )
ENDDO
! print THETA
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XTHT(IINDEX,JINDEX,JK) )
ENDDO
! print TKE
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XTKET(IINDEX,JINDEX,JK) )
ENDDO
! print RHODREF
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XRHODREF(IINDEX,JINDEX,JK) )
ENDDO
! print XDTHRAD
IF (CRAD .NE. "NONE") THEN
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XDTHRAD(IINDEX,JINDEX,JK) )
END DO
ELSE
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( 0.0 )
END DO
END IF
! print XCLDFR (the field is only allocated when radiation is used)
IF (CRAD .NE. "NONE") THEN
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XCLDFR(IINDEX,JINDEX,JK) )
END DO
ELSE DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( 0.0 )
END DO
END IF
! print JNO2
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
! the grid index of XJVALUES is starting with 1 for the physical grid
CALL WRITECLIP ( XJVALUES(IINDEX,JINDEX,JK-JPVEXT,2) )
ENDDO
! print JO1D
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
! the grid index of XJVALUES is starting with 1 for the physical grid
CALL WRITECLIP ( XJVALUES(IINDEX,JINDEX,JK-JPVEXT,3) )
ENDDO
! print water variables
DO JL = 1, NRR
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XRT(IINDEX,JINDEX,JK,JL) )
ENDDO
ENDDO
! print scalar variables
DO JL = 1, NSV
IF (JL>=NSV_CHEMBEG .AND. JL<=NSV_CHEMEND) THEN
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
! convert ppp to ppt
CALL WRITECLIP ( XSVT(IINDEX,JINDEX,JK,JL) * 1E12 )
ENDDO
ELSE
DO JK = NKMAX + JPVEXT , JPVEXT + 1, -1
CALL WRITECLIP ( XSVT(IINDEX,JINDEX,JK,JL) )
ENDDO
END IF
ENDDO
IF ((CPROGRAM =='DIAG ').AND.(LCHEMDIAG)) THEN
CALL IO_FILE_CLOSE_ll(TZFILE)
TZFILE => NULL()
END IF
IF (L1D) THEN
GSFIRSTCALL = .FALSE.
END IF
END IF
ENDDO
! increase timestep counter
ISCOUNT = ISCOUNT + 1
CONTAINS
SUBROUTINE WRITECLIP(PVAR)
USE MODD_CST, ONLY : XMNH_TINY
IMPLICIT NONE
REAL, INTENT(IN) :: PVAR
! clip very small numbers to zero, due to problems with PV-WAVE read
IF (ABS(PVAR) <= XMNH_TINY ) THEN
WRITE(ISIO1D,FMT=YFORM) 0.0
ELSE
WRITE(ISIO1D,FMT=YFORM) PVAR
END IF
RETURN
END SUBROUTINE WRITECLIP
END SUBROUTINE WRITE_TS1D