Newer
Older
!MNH_LIC Copyright 1994-2014 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.
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source$ $Revision$
! masdev4_7 BUG3 2007/11/20 12:25:23
!-----------------------------------------------------------------
! #################################
MODULE MODI_READ_ALL_DATA_GRIB_CASE
! #################################
INTERFACE
WAUTELET Philippe
committed
SUBROUTINE READ_ALL_DATA_GRIB_CASE(HFILE,TPPRE_REAL1,HGRIB,TPPGDFILE, &
WAUTELET Philippe
committed
USE MODD_IO_ll, ONLY: TFILEDATA
!
CHARACTER(LEN=4), INTENT(IN) :: HFILE ! which file ('ATM0','ATM1' or 'CHEM')
TYPE(TFILEDATA),POINTER,INTENT(INOUT) :: TPPRE_REAL1 ! PRE_REAL1 file
WAUTELET Philippe
committed
CHARACTER(LEN=28), INTENT(IN) :: HGRIB ! name of the GRIB file
TYPE(TFILEDATA), INTENT(IN) :: TPPGDFILE ! physiographic data file
INTEGER, INTENT(IN) :: KVERB ! verbosity level
LOGICAL, INTENT(IN) :: ODUMMY_REAL ! flag to interpolate dummy fields
REAL, INTENT(INOUT) :: PTIME_HORI ! time spent in hor. interpolations
!
END SUBROUTINE READ_ALL_DATA_GRIB_CASE
!
END INTERFACE
END MODULE MODI_READ_ALL_DATA_GRIB_CASE
! ##########################################################################
WAUTELET Philippe
committed
SUBROUTINE READ_ALL_DATA_GRIB_CASE(HFILE,TPPRE_REAL1,HGRIB,TPPGDFILE, &
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
PTIME_HORI,KVERB,ODUMMY_REAL )
! ##########################################################################
!
!!**** *READ_ALL_DATA_GRIB_CASE* - reads data for the initialization of real cases.
!!
!! PURPOSE
!! -------
! This routine reads the two input files :
! The PGD which is closed after reading
! The GRIB file
! Projection is read in READ_LFIFM_PGD (MODD_GRID).
! Grid and definition of large domain are read in PGD file and Grib files.
! The PGD files are also read in READ_LFIFM_PGD.
! The PGD file is closed.
! The MESO-NH domain is defined from PRE_REAL1.nam inputs in SET_SUBDOMAIN_CEP.
! Vertical grid is defined in READ_VER_GRID.
! PGD fields are stored on MESO-NH domain (in TRUNC_PGD).
!!
!!** METHOD
!! ------
!! 0. Declarations
!! 1. Declaration of arguments
!! 2. Declaration of local variables
!! 1. Read PGD file
!! 1. Domain restriction
!! 2. Coordinate conversion to lat,lon system
!! 2. Read Grib fields
!! 3. Vertical grid
!! 4. Free all temporary allocations
!!
!! EXTERNAL
!! --------
!! subroutine READ_LFIFM_PGD : to read PGD file
!! subroutine SET_SUBDOMAIN : to define the horizontal MESO-NH domain.
!! subroutine READ_VER_GRID : to read the vertical grid in namelist file.
!! subroutine HORIBL : horizontal bilinear interpolation
!! subroutine XYTOLATLON : projection from conformal to lat,lon
!!
!! Module MODI_SET_SUBDOMAIN : interface for subroutine SET_SUBDOMAIN
!! Module MODI_READ_VER_GRID : interface for subroutine READ_VER_GRID
!! Module MODI_HORIBL : interface for subroutine HORIBL
!! Module MODI_XYTOLATLON : interface for subroutine XYTOLATLON
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!!
!! Module MODD_CONF : contains configuration variables for all models.
!! NVERB : verbosity level for output-listing
!! Module MODD_LUNIT : contains logical unit names for all models
!! CLUOUT0 : name of output-listing
!! Module MODD_PGDDIM : contains dimension of PGD fields
!! NPGDIMAX: dimension along x (no external point)
!! NPGDJMAX: dimension along y (no external point)
!! Module MODD_PARAMETERS
!! JPHEXT
!!
!! REFERENCE
!! ---------
!!
!! Book 1 : Informations on ISBA model (soil moisture)
!! "Encoding and decoding Grib data", John D.Chambers, ECMWF, October 95
!! "A guide to Grib", John D.Stackpole, National weather service, March 94
!!
!! AUTHOR
!! ------
!!
!! J. Pettre and V. Bousquet
!!
!! MODIFICATIONS
!! -------------
!! Original 20/11/98
!! 15/03/99 (V. Masson) phasing with new PGD fields
!! 21/04/99 (V. Masson) bug in mask definitions for max Y index
!! 22/04/99 (V. Masson) optimizer bug in u,v loop
!! --> splitting of the loop
!! and splitting of the routine in more
!! contains
!! 28/05/99 (V. Bousquet) bug in wind interpolated variable for
!! Arpege
!! 31/05/99 (V. Masson) set pressure points (given on a regular grid at ECMWF)
!! on orography points (assuming the last are included in the former)
!! pressure computation from parameters A and B
!! (instead of interpolation from grib grid)
!! 20/07/00 (V. Masson) increase the threshold for land_sea index
!! 22/11/00 (P. Tulet) add INTERPOL_SV to initialize SV fields
!! (I. Mallet) from MOCAGE model (IMODE=3)
!! 01/02/01 (D. Gazen) add INI_NSV
!! 18/05/01 (P. Jabouille) problem with 129 grib code
!! 05/12/01 (I. Mallet) add Aladin reunion model
!! 02/10/02 (I. Mallet) 2 orography fields for CEP (SFC, ML=1)
!! 01/12/03 (D. Gazen) change Chemical scheme interface
!! 01/2004 (V. Masson) removes surface (externalization)
!! 01/06/02 (O.Nuissier) filtering of tropical cyclone
!! 01/05/04 (P. Tulet) add INTERPOL_SV to initialize SV dust
!! and aerosol fields
!! 08/06/2010 (G. Tanguy) replace GRIBEX by GRIB_API : change
!! of all the subroutine
!! 05/12/2016 (G.Delautier) length of HGRID for grib_api > 1.14
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
!------------
WAUTELET Philippe
committed
!
USE MODE_FM, ONLY: IO_FILE_CLOSE_ll
USE MODE_IO_ll, ONLY: UPCASE
WAUTELET Philippe
committed
USE MODE_THERMO
!
USE MODI_ADD_FORECAST_TO_DATE
USE MODI_READ_HGRID_n
USE MODI_READ_VER_GRID
USE MODI_XYTOLATLON
USE MODI_HORIBL
USE MODI_INI_NSV
USE MODI_CH_INIT_SCHEME_n
USE MODI_REMOVAL_VORTEX
USE MODI_CH_INIT_CCS
USE MODI_CH_AER_INIT_SOA
USE MODI_INI_CTURB
WAUTELET Philippe
committed
USE MODI_CH_OPEN_INPUT
WAUTELET Philippe
committed
USE MODD_IO_ll, ONLY: TFILEDATA
USE MODD_CONF
USE MODD_CONF_n
USE MODD_CST
USE MODD_LUNIT
USE MODD_PARAMETERS
USE MODD_GRID
USE MODD_GRID_n
USE MODD_DIM_n
USE MODD_PARAM_n, ONLY : CTURB
USE MODD_TIME
USE MODD_TIME_n
USE MODD_CH_MNHC_n, ONLY : LUSECHEM,LUSECHAQ,LUSECHIC,LCH_PH
USE MODD_CH_M9_n, ONLY : NEQ , CNAMES
USE MODD_CH_AEROSOL, ONLY: CORGANIC, NCARB, NSOA, NSP, LORILAM,&
JPMODE, LVARSIGI, LVARSIGJ
USE MODD_NSV , ONLY : NSV
USE MODD_HURR_CONF, ONLY : LFILTERING,CFILTERING
USE MODD_PREP_REAL
USE MODE_MODELN_HANDLER
!JUAN REALZ
USE MODE_MPPDB
!JUAN REALZ
WAUTELET Philippe
committed
!
USE GRIB_API
!
IMPLICIT NONE
!
!* 0.1. Declaration of arguments
! ------------------------
!
WAUTELET Philippe
committed
CHARACTER(LEN=4), INTENT(IN) :: HFILE ! which file ('ATM0','ATM1' or 'CHEM')
TYPE(TFILEDATA),POINTER,INTENT(INOUT) :: TPPRE_REAL1! PRE_REAL1 file
WAUTELET Philippe
committed
CHARACTER(LEN=28), INTENT(IN) :: HGRIB ! name of the GRIB file
TYPE(TFILEDATA), INTENT(IN) :: TPPGDFILE ! physiographic data file
INTEGER, INTENT(IN) :: KVERB ! verbosity level
LOGICAL, INTENT(IN) :: ODUMMY_REAL ! flag to interpolate dummy fields
REAL, INTENT(INOUT) :: PTIME_HORI ! time spent in hor. interpolations
!
!* 0.2 Declaration of local variables
! ------------------------------
! General purpose variables
INTEGER :: ILUOUT0 ! Unit used for output msg.
INTEGER :: IRESP ! Return code of FM-routines
INTEGER :: IRET ! Return code from subroutines
WAUTELET Philippe
committed
INTEGER(KIND=kindOfInt) :: IRET_GRIB ! Return code from subroutines
REAL :: ZA,ZB,ZC ! Dummy variables
REAL :: ZD,ZE,ZF ! |
REAL :: ZTEMP ! |
INTEGER :: JI,JJ ! Dummy counters
INTEGER :: JLOOP1,JLOOP2 ! |
INTEGER :: JLOOP3,JLOOP4 ! |
INTEGER :: JLOOP ! |
! Variables used by the PGD reader
CHARACTER(LEN=28) :: YPGD_NAME ! not used - dummy argument
CHARACTER(LEN=28) :: YPGD_DAD_NAME ! not used - dummy argument
CHARACTER(LEN=2) :: YPGD_TYPE ! not used - dummy argument
! PGD Grib definition variables
INTEGER :: INO ! Number of points of the grid
INTEGER :: IIU ! Number of points along X
INTEGER :: IJU ! Number of points along Y
REAL, DIMENSION(:), ALLOCATABLE :: ZXOUT ! mapping PGD -> Grib (lon.)
REAL, DIMENSION(:), ALLOCATABLE :: ZYOUT ! mapping PGD -> Grib (lat.)
REAL, DIMENSION(:), ALLOCATABLE :: ZLONOUT ! mapping PGD -> Grib (lon.)
REAL, DIMENSION(:), ALLOCATABLE :: ZLATOUT ! mapping PGD -> Grib (lat.)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZXM ! X of PGD mass points
REAL, DIMENSION(:,:), ALLOCATABLE :: ZYM ! Y of PGD mass points
REAL, DIMENSION(:,:), ALLOCATABLE :: ZLATM ! Lat of PGD mass points
REAL, DIMENSION(:,:), ALLOCATABLE :: ZLONM ! Lon of PGD mass points
! Variable involved in the task of reading the grib file
INTEGER(KIND=kindOfInt) :: IUNIT ! unit of the grib file
CHARACTER(LEN=50) :: HGRID ! type of grid
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
INTEGER :: IPAR ! Parameter identifier
INTEGER :: ITYP ! type of level (Grib code table 3)
INTEGER :: ILEV1 ! level definition
INTEGER :: ILEV2 ! level definition
INTEGER :: IMODEL ! Type of Grib file :
! 0 -> ECMWF
! 1 -> METEO FRANCE - ALADIN/AROME
! 2 -> METEO FRANCE - ALADIN-REUNION
! 3 -> METEO FRANCE - ARPEGE
! 4 -> METEO FRANCE - ARPEGE
! 5 -> METEO FRANCE - MOCAGE
INTEGER :: ICENTER ! number of center
INTEGER :: ISIZE ! size of grib message
INTEGER(KIND=kindOfInt) :: ICOUNT ! number of messages in the file
INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE :: IGRIB ! number of the grib in memory
INTEGER :: INUM ,INUM_ZS ! number of a grib message
REAL,DIMENSION(:),ALLOCATABLE :: ZPARAM ! parameter of girb grid
INTEGER,DIMENSION(:),ALLOCATABLE :: IINLO ! longitude of grib grid
INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE :: IINLO_GRIB ! longitude of grib grid
REAL,DIMENSION(:),ALLOCATABLE :: ZPARAM_ZS ! parameter of girb grid for ZS
INTEGER,DIMENSION(:),ALLOCATABLE :: IINLO_ZS ! longitude of grib grid for ZS
REAL,DIMENSION(:),ALLOCATABLE :: ZVALUE ! Intermediate array
REAL,DIMENSION(:),ALLOCATABLE :: ZOUT ! Intermediate arrays
! Grib grid definition variables
INTEGER :: INI ! Number of points
INTEGER :: INLEVEL ! Number of levels
INTEGER :: ISTARTLEVEL ! First level (0 or 1)
TYPE(DATE_TIME) :: TPTCUR ! Date & time of the grib file data
INTEGER :: ITWOZS
! surface pressure
REAL, DIMENSION(:), ALLOCATABLE :: ZPS_G ! Grib data : Ps
REAL, DIMENSION(:), ALLOCATABLE :: ZLNPS_G ! Grib data : ln(Ps)
REAL, DIMENSION(:), ALLOCATABLE :: ZWORK_LNPS ! Grib data on zs grid: ln(Ps)
INTEGER :: INJ,INJ_ZS
! orography
CHARACTER(LEN=50) :: HGRID_ZS ! type of grid
!
! Reading and projection of the wind vectors u, v
REAL :: ZALPHA ! Angle of rotation
REAL, DIMENSION(:), ALLOCATABLE :: ZTU_LS ! Intermediate array for U
REAL, DIMENSION(:), ALLOCATABLE :: ZTV_LS ! | V
REAL :: ZLATPOLE ! Arpege stretching pole latitude
REAL :: ZLONPOLE ! Arpege stretching pole longitude
REAL :: ZLAT,ZLON ! Lat,lon of current point
REAL :: ZCOS,ZSIN ! cos,sin of rotation matrix
REAL, DIMENSION(:), ALLOCATABLE :: ZTU0_LS ! Arpege temp array for U
REAL, DIMENSION(:), ALLOCATABLE :: ZTV0_LS ! | V
!
! variables for hurricane filtering
REAL, DIMENSION(:,:), ALLOCATABLE :: ZTVF_LS,ZMSLP_LS
REAL :: ZGAMREF ! Standard atmosphere lapse rate (K/m)
! date
INTEGER :: ITIME
INTEGER :: IDATE
INTEGER :: ITIMESTEP
CHARACTER(LEN=10) :: CSTEPUNIT
!chemistery field
WAUTELET Philippe
committed
CHARACTER(LEN=16) :: YPRE_MOC="PRE_MOC1.nam"
INTEGER, DIMENSION(:), ALLOCATABLE :: INUMGRIB, INUMLEV ! grib
INTEGER, DIMENSION(:), ALLOCATABLE :: INUMLEV1, INUMLEV2 !numbers
INTEGER :: IMOC
INTEGER :: IVAR
INTEGER :: ICHANNEL
INTEGER :: INDX
INTEGER :: INACT
CHARACTER(LEN=40) :: YINPLINE ! input line
WAUTELET Philippe
committed
CHARACTER(LEN=16) :: YFIELD
CHARACTER, PARAMETER :: YPTAB = CHAR(9) ! TAB character is ASCII : 9
CHARACTER, PARAMETER :: YPCOM = CHAR(44)! COMma character is ASCII : 44
CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE :: YMNHNAME ! species names
INTEGER :: JN, JNREAL ! loop control variables
CHARACTER(LEN=40) :: YFORMAT
! temperature and humidity
INTEGER :: IT,IQ
REAL, DIMENSION(:,:), ALLOCATABLE :: ZPF_G ! Pressure (flux point)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZPM_G ! Pressure (mass point)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZEXNF_G ! Exner fct. (flux point)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZEXNM_G ! Exner fct. (mass point)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZT_G ! Temperature
REAL, DIMENSION(:,:), ALLOCATABLE :: ZQ_G ! Specific humidity
REAL, DIMENSION(:), ALLOCATABLE :: ZH_G ! Relative humidity
REAL, DIMENSION(:), ALLOCATABLE :: ZTHV_G ! Theta V
REAL, DIMENSION(:), ALLOCATABLE :: ZRV_G ! Vapor mixing ratio
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZPF_LS ! Pressure (flux point)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZPM_LS ! Pressure (mass point)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEXNF_LS ! Exner fct. (flux point)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEXNM_LS ! Exner fct. (mass point)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZH_LS ! Relative humidity
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRV_LS ! Vapor mixing ratio
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTHV_LS ! Theta V
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTEV_LS ! T V
REAL, DIMENSION(:), ALLOCATABLE :: ZPV ! vertical level in grib file
INTEGER :: IPVPRESENT ,IPV
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZR_DUM
INTEGER :: IMI
WAUTELET Philippe
committed
TYPE(TFILEDATA),POINTER :: TZFILE
!
!---------------------------------------------------------------------------------------
!
WAUTELET Philippe
committed
TZFILE => NULL()
!
IMI = GET_CURRENT_MODEL_INDEX()
!
!* 1. READ PGD FILE
! -------------
!
WAUTELET Philippe
committed
CALL READ_HGRID_n(TPPGDFILE,YPGD_NAME,YPGD_DAD_NAME,YPGD_TYPE)
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
!
! 1.1 Domain restriction
!
!JUAN REALZ
CALL GET_DIM_EXT_ll('B',IIU,IJU)
!!$IIU=NIMAX + 2*JPHEXT
!!$IJU=NJMAX + 2*JPHEXT
!JUAN REALZ
INO = IIU * IJU
!
!
! 1.2 Coordinate conversion to lat,lon system
!
ALLOCATE (ZXM(IIU,IJU))
ALLOCATE (ZYM(IIU,IJU))
ALLOCATE (ZLONM(IIU,IJU))
ALLOCATE (ZLATM(IIU,IJU))
ZXM(1:IIU-1,1) = (XXHAT(1:IIU-1) + XXHAT(2:IIU) ) / 2.
ZXM(IIU,1) = XXHAT(IIU) - XXHAT(IIU-1) + ZXM(IIU-1,1)
ZXM(:,2:IJU) = SPREAD(ZXM(:,1),2,IJU-1)
ZYM(1,1:IJU-1) = (XYHAT(1:IJU-1) + XYHAT(2:IJU)) / 2.
ZYM(1,IJU) = XYHAT(IJU) - XYHAT(IJU-1) + ZYM(1,IJU-1)
ZYM(2:IIU,:) = SPREAD(ZYM(1,:),1,IIU-1)
CALL SM_XYTOLATLON_A (XLAT0,XLON0,XRPK,XLATORI,XLONORI,ZXM,ZYM,ZLATM,ZLONM, &
IIU,IJU)
ALLOCATE (ZLONOUT(INO))
ALLOCATE (ZLATOUT(INO))
JLOOP1 = 0
DO JJ = 1, IJU
ZLONOUT(JLOOP1+1:JLOOP1+IIU) = ZLONM(1:IIU,JJ)
ZLATOUT(JLOOP1+1:JLOOP1+IIU) = ZLATM(1:IIU,JJ)
JLOOP1 = JLOOP1 + IIU
ENDDO
DEALLOCATE (ZLATM)
DEALLOCATE (ZLONM)
DEALLOCATE (ZYM)
DEALLOCATE (ZXM)
!
ALLOCATE (ZXOUT(INO))
ALLOCATE (ZYOUT(INO))
!
!---------------------------------------------------------------------------------------
!
!* 2. READ GRIB FIELDS
! ----------------
!
IF (HFILE(1:3)=='ATM' .OR. HFILE=='CHEM') THEN
WRITE (ILUOUT0,'(A,A4)') ' -- Grib reader started for ',HFILE
ELSE
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','bad input argument')
END IF
!
!* 2.1 Charge in memory the grib messages
!
! open grib file
CALL GRIB_OPEN_FILE(IUNIT,HGRIB,'R',IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error opening the grib file ',TRIM(HGRIB),', error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
! count the messages in the file
CALL GRIB_COUNT_IN_FILE(IUNIT,ICOUNT,IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error in reading the grib file ',TRIM(HGRIB),', error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
ALLOCATE(IGRIB(ICOUNT))
! initialize the tabular with a negativ number
! ( all the IGRIB will be different )
IGRIB(:)=-12
!charge all the message in memory
DO JLOOP=1,ICOUNT
CALL GRIB_NEW_FROM_FILE(IUNIT,IGRIB(JLOOP),IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error in reading the grib file - ILOOP=',JLOOP,' - error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
END DO
! close the grib file
CALL GRIB_CLOSE_FILE(IUNIT)
!
!
!---------------------------------------------------------------------------------------
!* 2.2 Research center with the first message
!---------------------------------------------------------------------------------------
!
CALL GRIB_GET(IGRIB(1),'centre',ICENTER,IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error in reading center - error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
CALL GRIB_GET(IGRIB(1),'typeOfGrid',HGRID,IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error in reading type of grid - error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
!
IMODEL = -1
SELECT CASE (ICENTER)
CASE (98)
WRITE (ILUOUT0,'(A)') &
' | Grib file from European Center for Medium-range Weather Forecast'
IMODEL = 0
ALLOCATE(ZPARAM(6))
CASE (85)
SELECT CASE (HGRID)
CASE('lambert')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Aladin france model'
IMODEL = 1
ALLOCATE(ZPARAM(10))
CASE('mercator')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Aladin reunion model'
IMODEL = 2
ALLOCATE(ZPARAM(9))
CASE('unknown_PLPresent','reduced_stretched_rotated_gg')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arpege model'
IMODEL = 3
ALLOCATE(ZPARAM(10))
CASE('regular_gg')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arpege model'
WRITE (ILUOUT0,'(A)') 'but same grid as ECMWF model (unstretched)'
IMODEL = 4
ALLOCATE(ZPARAM(10))
CASE('regular_ll')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Mocage model'
IMODEL = 5
ALLOCATE(ZPARAM(6))
END SELECT
END SELECT
IF (IMODEL==-1) THEN
!callabortstop
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','unsupported Grib file format')
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
END IF
!
!---------------------------------------------------------------------------------------
!* 2.3 Read and interpol orography
!---------------------------------------------------------------------------------------
!
WRITE (ILUOUT0,'(A)') ' | Searching orography'
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
CALL SEARCH_FIELD(129,109,1,0,IGRIB,INUM_ZS)
IF(INUM_ZS < 0) THEN
CALL SEARCH_FIELD(129,1,0,0,IGRIB,INUM_ZS)
IF(INUM_ZS < 0) THEN
WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
END IF
END IF
CASE(3,4,5) ! arpege et mocage
CALL SEARCH_FIELD(8,-1,-1,-1,IGRIB,INUM_ZS)
IF(INUM_ZS < 0) THEN
WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
ENDIF
CASE(1,2) ! aladin et aladin reunion
CALL SEARCH_FIELD(6,-1,-1,-1,IGRIB,INUM_ZS)
IF(INUM_ZS < 0) THEN
WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
ENDIF
END SELECT
ZPARAM(:)=-999.
CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ,IRET_GRIB)
ALLOCATE(IINLO(INJ))
CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM_ZS),IIU,IJU,ZLONOUT,ZLATOUT,&
ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
ALLOCATE(ZPARAM_ZS(SIZE(ZPARAM)))
ZPARAM_ZS=ZPARAM
ALLOCATE(IINLO_ZS(SIZE(IINLO)))
IINLO_ZS=IINLO
CALL GRIB_GET_SIZE(IGRIB(INUM_ZS),'values',ISIZE)
ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM_ZS),'values',ZVALUE)
ALLOCATE(ZOUT(INO))
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
DEALLOCATE(IINLO)
DEALLOCATE(ZVALUE)
!
! Datas given in archives are multiplied by the gravity acceleration
ZOUT = ZOUT / XG
!
! Stores the field in a 2 dimension array
IF (HFILE(1:3)=='ATM') THEN
ALLOCATE (XZS_LS(IIU,IJU))
ALLOCATE (XZSMT_LS(IIU,IJU))
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XZS_LS)
XZSMT_LS = XZS_LS ! no smooth orography in this case
ELSE IF (HFILE=='CHEM') THEN
ALLOCATE (XZS_SV_LS(IIU,IJU))
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XZS_SV_LS)
END IF
DEALLOCATE (ZOUT)
!
!---------------------------------------------------------------------------------------
!* 2.4 Interpolation surface pressure
!---------------------------------------------------------------------------------------
!
!* 2.4.1 Read pressure
!
WRITE (ILUOUT0,'(A)') ' | Searching pressure'
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
CALL SEARCH_FIELD(152,-1,-1,-1,IGRIB,INUM)
CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
CALL SEARCH_FIELD(1,-1,-1,-1,IGRIB,INUM)
END SELECT
IF(INUM < 0) THEN
WRITE (ILUOUT0,'(A)')'Surface pressure is missing - abort'
CALL ABORT
STOP
ENDIF
! recuperation du tableau de valeurs
CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
! determination des tableaux ZPS_G et ZLNPS_G
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
ALLOCATE (ZPS_G (ISIZE))
ALLOCATE (ZLNPS_G(ISIZE))
ZLNPS_G(:) = ZVALUE(1:ISIZE)
ZPS_G (:) = EXP(ZVALUE(1:ISIZE))
CASE(1,2,3,4,5) ! arpege mocage aladin aladin-reunion
ALLOCATE (ZPS_G (ISIZE))
ALLOCATE (ZLNPS_G(ISIZE))
ZPS_G (:) = ZVALUE(1:ISIZE)
ZLNPS_G(:) = LOG(ZVALUE(1:ISIZE))
END SELECT
DEALLOCATE (ZVALUE)
!
!* 2.4.2 Removes pressure points not on orography points
! (if pressure is on a regular grid)
CALL GRIB_GET(IGRIB(INUM),'typeOfGrid',HGRID)
CALL GRIB_GET(IGRIB(INUM_ZS),'typeOfGrid',HGRID_ZS)
CALL GRIB_GET(IGRIB(INUM),'Nj',INJ)
CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ_ZS)
!
IF ( HGRID(1:7)=='regular' .AND. HGRID_ZS(1:7)=='reduced' .AND.&
INJ == INJ_ZS) THEN
WRITE (ILUOUT0,'(A)')'HGRID(1:7)==regular .AND. HGRID_ZS(1:7)==reduced .AND. INJ == INJ_ZS - abort'
CALL ABORT
STOP
ELSE
ALLOCATE(ZWORK_LNPS(SIZE(ZLNPS_G)))
ZWORK_LNPS(:) = ZLNPS_G(:)
ENDIF
!
IF (HFILE(1:3)=='ATM') THEN
ALLOCATE (XPS_LS(IIU,IJU))
ELSE IF (HFILE=='CHEM') THEN
ALLOCATE (XPS_SV_LS(IIU,IJU))
END IF
!
ALLOCATE(IINLO(INJ))
CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
ALLOCATE(ZOUT(INO))
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI,&
ZWORK_LNPS,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE. )
DEALLOCATE(ZWORK_LNPS)
DEALLOCATE(IINLO)
!
!* 2.4.3 conversion to surface pressure
!
ZOUT=EXP(ZOUT)
IF (HFILE(1:3)=='ATM') THEN
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XPS_LS(:,:))
ELSE IF (HFILE=='CHEM') THEN
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XPS_SV_LS(:,:))
END IF
!JUAN REALZ
CALL MPPDB_CHECK2D(XZS_LS,"XZS_LS",PRECISION)
!JUAN REALZ
DEALLOCATE (ZOUT)
DEALLOCATE (ZLNPS_G)
!
!---------------------------------------------------------------------------------------
!* 2.5 Interpolation temperature and humidity
!---------------------------------------------------------------------------------------
!
!* 2.5.1 Read T and Q on each level
!
WRITE (ILUOUT0,'(A)') ' | Reading T and Q fields'
!
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
ISTARTLEVEL=1
IT=130
IQ=133
CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
IF(INUM < 0) THEN
ISTARTLEVEL=0
CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
ENDIF
IF(INUM < 0) THEN
WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
CALL ABORT
STOP
ENDIF
CALL SEARCH_FIELD(IQ,109,ISTARTLEVEL,-1,IGRIB,INUM)
IF(INUM < 0) THEN
WRITE (ILUOUT0,'(A)')'Atmospheric specific humidity is missing - abort'
CALL ABORT
STOP
ENDIF
CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
IT=11
IQ=51
ISTARTLEVEL=1
CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
IF(INUM < 0) THEN
ISTARTLEVEL=0
CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
ENDIF
IF(INUM < 0) THEN
WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
CALL ABORT
STOP
ENDIF
CALL SEARCH_FIELD(IQ,109,ISTARTLEVEL,-1,IGRIB,INUM)
IF(INUM < 0) THEN
WRITE (ILUOUT0,'(A)')'Atmospheric specific humidity is missing - abort'
CALL ABORT
STOP
ENDIF
END SELECT
!
CALL GRIB_GET(IGRIB(INUM),'NV',INLEVEL)
INLEVEL = NINT(INLEVEL / 2.) - 1
CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
!
ALLOCATE (ZT_G(ISIZE,INLEVEL))
ALLOCATE (ZQ_G(ISIZE,INLEVEL))
!
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1-1+ISTARTLEVEL
CALL SEARCH_FIELD(IQ,109,ILEV1,-1,IGRIB,INUM)
IF (INUM< 0) THEN
!callabortstop
WRITE(YMSG,*) 'atmospheric humidity level ',JLOOP1,' is missing'
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
CALL GRIB_GET(IGRIB(INUM),'values',ZQ_G(:,INLEVEL-JLOOP1+1))
CALL SEARCH_FIELD(IT,109,ILEV1,-1,IGRIB,INUM)
IF (INUM< 0) THEN
!callabortstop
WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
END IF
CALL GRIB_GET(IGRIB(INUM),'values',ZT_G(:,INLEVEL-JLOOP1+1))
CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
END DO
ALLOCATE(IINLO(INJ))
CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
!
!* 2.5.2 Load level definition parameters A and B
!
IF (HFILE(1:3)=='ATM') THEN
XP00_LS = 101325.
ELSE IF (HFILE=='CHEM') THEN
XP00_SV_LS = 101325.
END IF
!
IF (INLEVEL > 0) THEN
IF (HFILE(1:3)=='ATM') THEN
ALLOCATE (XA_LS(INLEVEL))
ALLOCATE (XB_LS(INLEVEL))
ELSE IF (HFILE=='CHEM') THEN
ALLOCATE (XA_SV_LS(INLEVEL))
ALLOCATE (XB_SV_LS(INLEVEL))
END IF
!
CALL GRIB_GET(IGRIB(INUM),'PVPresent',IPVPRESENT)
IF (IPVPRESENT==1) THEN
CALL GRIB_GET_SIZE(IGRIB(INUM),'pv',IPV)
ALLOCATE(ZPV(IPV))
CALL GRIB_GET(IGRIB(INUM),'pv',ZPV)
ELSE
!callabortstop
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','there is no PV value in this message')
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
ENDIF
SELECT CASE (IMODEL)
CASE (0,3,4)
DO JLOOP1 = 1, INLEVEL
XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+JLOOP1) / XP00_LS
XB_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
END DO
CASE (1,2)
JLOOP2 = 2
DO JLOOP1 = 1, INLEVEL
JLOOP2 = JLOOP2 + 1
XA_LS(1 + INLEVEL - JLOOP1) = ZPV(JLOOP2)
JLOOP2 = JLOOP2 + 1
XB_LS(1 + INLEVEL - JLOOP1) = ZPV(JLOOP2)
END DO
CASE (5)
DO JLOOP1 = 1, INLEVEL
IF (HFILE(1:3)=='ATM') THEN
XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+ JLOOP1) / XP00_LS**2
XB_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
ELSE IF (HFILE=='CHEM') THEN
XA_SV_LS(1 + INLEVEL - JLOOP1) = ZPV(1+ JLOOP1) / XP00_LS**2
XB_SV_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
END IF
END DO
END SELECT
ELSE
!callabortstop
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','level definition section is missing')
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
END IF
!
!* 2.5.3 Compute atmospheric pressure on grib grid
!
WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on Grib grid is being computed'
ALLOCATE (ZPF_G(INI,INLEVEL))
IF (HFILE(1:3)=='ATM') THEN
ZPF_G(:,:) = SPREAD(XA_LS,1,INI)*XP00_LS + &
SPREAD(XB_LS,1,INI)*SPREAD(ZPS_G(1:INI),2,INLEVEL)
ELSE IF (HFILE=='CHEM') THEN
ZPF_G(:,:) = SPREAD(XA_SV_LS,1,INI)*XP00_SV_LS + &
SPREAD(XB_SV_LS,1,INI)*SPREAD(ZPS_G(1:INI),2,INLEVEL)
END IF
DEALLOCATE (ZPS_G)
!
ALLOCATE (ZEXNF_G(INI,INLEVEL))
ZEXNF_G(:,:) = (ZPF_G(:,:)/XP00)**(XRD/XCPD)
!
ALLOCATE (ZEXNM_G(INI,INLEVEL))
ZEXNM_G(:,1:INLEVEL-1) = (ZEXNF_G(:,1:INLEVEL-1)-ZEXNF_G(:,2:INLEVEL)) / &
(LOG(ZEXNF_G(:,1:INLEVEL-1))-LOG(ZEXNF_G(:,2:INLEVEL)))
ZEXNM_G(:,INLEVEL) = (ZPF_G(:,INLEVEL)/2./XP00)**(XRD/XCPD)
DEALLOCATE (ZEXNF_G)
DEALLOCATE (ZPF_G)
!
ALLOCATE (ZPM_G(INI,INLEVEL))
ZPM_G(:,:) = XP00*(ZEXNM_G(:,:))**(XCPD/XRD)
!
!* 2.5.4 Compute the relative humidity and the interpolate Thetav, P, Q and H
!
IF (IMODEL==1) THEN
! search cloud_water in Arome case (forecast)
ISTARTLEVEL = 1
IPAR=246
CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
IF (INUM < 0) THEN
ISTARTLEVEL = 0
CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
END IF
IF (INUM > 0) THEN
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arome model (forecast)'
LCPL_AROME=.TRUE.
NRR=6
END IF
! search also turbulent kinetic energy
ISTARTLEVEL = 1
IPAR=251
CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
IF (INUM < 0) THEN
ISTARTLEVEL = 0
CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
END IF
IF (INUM > 0) CTURB='TKEL'
END IF
!
!
WRITE (ILUOUT0,'(A)') ' | Computing relative humidity on each level'
ALLOCATE (ZH_G(INI))
ALLOCATE (ZH_LS(IIU,IJU,INLEVEL))
IF (HFILE(1:3)=='ATM') THEN
ALLOCATE (XT_LS(IIU,IJU,INLEVEL))
ALLOCATE (XQ_LS(IIU,IJU,INLEVEL,NRR)) ; XQ_LS=0.
ELSE IF (HFILE=='CHEM') THEN
ALLOCATE (XT_SV_LS(IIU,IJU,INLEVEL))
ALLOCATE (XQ_SV_LS(IIU,IJU,INLEVEL,1))
END IF
IF (CTURB=='TKEL') THEN
IF (ALLOCATED(XTKE_LS)) DEALLOCATE(XTKE_LS)
ALLOCATE(XTKE_LS(IIU,IJU,INLEVEL)) ; XTKE_LS=0.
CALL INI_CTURB
ELSE
IF (ALLOCATED(XTKE_LS)) DEALLOCATE(XTKE_LS)
ALLOCATE(XTKE_LS(0,0,0))
END IF
ALLOCATE (ZTHV_LS(IIU,IJU,INLEVEL))
ALLOCATE (ZTHV_G(INI))
ALLOCATE (ZRV_G(INI))
ALLOCATE (ZOUT(INO))
DO JLOOP1=1, INLEVEL
!
! Compute Theta V and relative humidity on grib grid
!
! (1/rv) = (1/q) - 1
! Thetav = T . (P0/P)^(Rd/Cpd) . ( (1 + (Rv/Rd).rv) / (1 + rv) )
! Hu = P / ( ( (Rd/Rv) . ((1/rv) - 1) + 1) . Es(T) )
!
ZRV_G(:) = 1. / (1./MAX(ZQ_G(:,JLOOP1),1.E-12) - 1.)
!
ZTHV_G(:)=ZT_G(:,JLOOP1) * ((XP00/ZPM_G(:,JLOOP1))**(XRD/XCPD)) * &
((1. + XRV*ZRV_G(:)/XRD) / (1. + ZRV_G(:)))
!
ZH_G(1:INI) = 100. * ZPM_G(:,JLOOP1) / ( (XRD/XRV)*(1./ZRV_G(:)+1.)*SM_FOES(ZT_G(:,JLOOP1)) )
ZH_G(:) = MAX(MIN(ZH_G(:),100.),0.)
!
! Interpolation : H
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
ZH_G,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZH_LS(:,:,JLOOP1))
ZH_LS(:,:,JLOOP1) = MAX(MIN(ZH_LS(:,:,JLOOP1),100.),0.)
!
! interpolation : Theta V
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
ZTHV_G,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZTHV_LS(:,:,JLOOP1))
!
END DO
DEALLOCATE (ZOUT)
!
!* 2.5.5 Compute atmospheric pressure on MESO-NH grid
!
WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on MesoNH grid is being computed'
ALLOCATE (ZPF_LS(IIU,IJU,INLEVEL))
IF (HFILE(1:3)=='ATM') THEN
ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_LS,1,IIU),2,IJU)*XP00_LS + &
SPREAD(SPREAD(XB_LS,1,IIU),2,IJU)*SPREAD(XPS_LS,3,INLEVEL)
ELSE IF (HFILE=='CHEM') THEN
ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_SV_LS,1,IIU),2,IJU)*XP00_LS + &
SPREAD(SPREAD(XB_SV_LS,1,IIU),2,IJU)*SPREAD(XPS_SV_LS,3,INLEVEL)
END IF
!
ALLOCATE (ZEXNF_LS(IIU,IJU,INLEVEL))
ZEXNF_LS(:,:,:) = (ZPF_LS(:,:,:)/XP00)**(XRD/XCPD)
!
ALLOCATE (ZEXNM_LS(IIU,IJU,INLEVEL))
ZEXNM_LS(:,:,1:INLEVEL-1) = (ZEXNF_LS(:,:,1:INLEVEL-1)-ZEXNF_LS(:,:,2:INLEVEL)) / &
(LOG(ZEXNF_LS(:,:,1:INLEVEL-1))-LOG(ZEXNF_LS(:,:,2:INLEVEL)))
ZEXNM_LS(:,:,INLEVEL) = (ZPF_LS(:,:,INLEVEL)/2./XP00)**(XRD/XCPD)
DEALLOCATE (ZEXNF_LS)
DEALLOCATE (ZPF_LS)
!
ALLOCATE (ZPM_LS(IIU,IJU,INLEVEL))
ZPM_LS(:,:,:) = XP00*(ZEXNM_LS(:,:,:))**(XCPD/XRD)
!
!* 2.5.6 Compute the vapor mixing ratio and the final specific humdity
!
! The vapor mixing ratio is calculated by an interating process on rv and
! Thetav. Have a look to MODE_THERMO for further informations.
ALLOCATE (ZR_DUM(IIU,IJU,INLEVEL,1))
ALLOCATE (ZRV_LS(IIU,IJU,INLEVEL))
ALLOCATE (ZTEV_LS(IIU,IJU,INLEVEL))
ZTEV_LS(:,:,:) = ZTHV_LS(:,:,:) * ZEXNM_LS(:,:,:)
ZRV_LS(:,:,:) = SM_PMR_HU(ZPM_LS(:,:,:), &
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
ZTEV_LS(:,:,:),ZH_LS(:,:,:),ZR_DUM(:,:,:,:),KITERMAX=100)
IF (HFILE(1:3)=='ATM') THEN
XQ_LS(:,:,:,1) = ZRV_LS(:,:,:) / (1. + ZRV_LS(:,:,:))
ELSE IF (HFILE=='CHEM') THEN
XQ_SV_LS(:,:,:,1) = ZRV_LS(:,:,:) / (1. + ZRV_LS(:,:,:))
ENDIF
!JUAN
CALL MPPDB_CHECK3D(XQ_LS(:,:,:,1),"XQ_LS",PRECISION)
!JUAN
DEALLOCATE (ZTEV_LS)
DEALLOCATE (ZH_LS)
DEALLOCATE (ZR_DUM)
!
!* 2.5.7 Compute T from the interpolated Theta V
!
! T = Thetav . (P/P0)^(Rd/Cpd) . ((1 + rv) / (1 + (Rv/Rd).rv))
!!
IF (HFILE(1:3)=='ATM') THEN
XT_LS(:,:,:) = ZTHV_LS(:,:,:) * ZEXNM_LS(:,:,:) &
* ((1.+ZRV_LS(:,:,:))/(1.+(XRV/XRD)*ZRV_LS(:,:,:)))
CALL MPPDB_CHECK3D(XT_LS,"XT_LS",PRECISION)
ELSE IF (HFILE=='CHEM') THEN
XT_SV_LS(:,:,:) = ZTHV_LS(:,:,:) * ZEXNM_LS(:,:,:) &
* ((1.+ZRV_LS(:,:,:))/(1.+(XRV/XRD)*ZRV_LS(:,:,:)))
CALL MPPDB_CHECK3D(XT_SV_LS,"XT_SV_LS",PRECISION)
ENDIF
!
DEALLOCATE (ZRV_LS)
DEALLOCATE (ZTHV_LS)
DEALLOCATE (ZPM_LS)
DEALLOCATE (ZEXNM_LS)
!
!* 2.5.8 Read the other specific ratios (if Arome model)
!
IF (NRR >1) THEN
WRITE (ILUOUT0,'(A)') ' | Reading Q fields (except humidity)'
DO JLOOP2=1,NRR-1
IPAR=246+JLOOP2-1
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1-1+ISTARTLEVEL
CALL SEARCH_FIELD(IPAR,109,ILEV1,-1,IGRIB,INUM)
IF (INUM < 0) THEN
!callabortstop
WRITE(YMSG,*) 'Specific ratio ',IPAR,' at level ',JLOOP1,' is missing'
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
ALLOCATE(ZOUT(INO))
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,1+JLOOP2))
DEALLOCATE (ZVALUE)
DEALLOCATE (ZOUT)
END DO
END DO
END IF
!
IF (CTURB=='TKEL') THEN
WRITE (ILUOUT0,'(A)') ' | Reading TKE field'
IPAR=251
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1-1+ISTARTLEVEL
CALL SEARCH_FIELD(IPAR,109,ILEV1,-1,IGRIB,INUM)
IF (INUM < 0) THEN
!callabortstop
WRITE(YMSG,*) 'TKE at level ',JLOOP1,' is missing'
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
END IF
CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
ALLOCATE(ZOUT(INO))
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XTKE_LS(:,:,INLEVEL-JLOOP1+1))
DEALLOCATE (ZVALUE)
DEALLOCATE (ZOUT)
END DO
END IF
DEALLOCATE(IINLO)
!
!---------------------------------------------------------------------------------------
!* 2.6 Interpolation of MOCAGE variable
!---------------------------------------------------------------------------------------
IF (IMODEL==5) THEN
! Always initialize chemical scheme variables before INI_NSV call !
CALL CH_INIT_SCHEME_n(IMI,LUSECHAQ,LUSECHIC,LCH_PH,ILUOUT0,KVERB)
LUSECHEM = .TRUE.
IF (LORILAM) THEN
CORGANIC = "MPMPO"
LVARSIGI = .TRUE.
LVARSIGJ = .TRUE.
CALL CH_AER_INIT_SOA(ILUOUT0, KVERB)
END IF
! initialise NSV_* variables
CALL INI_NSV(1)
IF( HFILE=='ATM0' ) THEN
ALLOCATE (XSV_LS(IIU,IJU,INLEVEL,NSV))
ELSE IF (HFILE=='CHEM' ) THEN
DEALLOCATE(XSV_LS)
ALLOCATE (XSV_LS(IIU,IJU,INLEVEL,NSV))