Ignore:
Timestamp:
Jul 31, 2015, 7:22:21 PM (9 years ago)
Author:
dcugnet
Message:
  • Add parallel capability for ce0l.
  • Small bug in grid_noro fixed (smoothed topography was used instead of unsmoothed one for geopotential computation at north pole).
  • Removed average of mass at poles in etat0dyn_netcdf after start_init_dyn => different results in the zoomed grid case.
  • ok_etat0=n and ok_limit=y combination now works fine (if no initial state is needed, but only limit.nc file). This required:
    • to move grid_noro0 and start_init_noro0 subroutines from etat0dyn_netcdf.F90 to limit_netcdf.F90
    • to create init_ssrf_m.F90 file, so that sub-surfaces can be initialized from limit_netcdf.F90 without any etat0*_netcdf routines call).
  • Simplified somehow the corresponding code, in particular: 1) removed obsolete flags "oldice". 2) removed flag "ibar": barycentric interpolation is used everywhere (except in start_init_subsurf, still calling grille_m - to be changed soon). 3) removed useless CPP_PHY precompilation directives, considering the possibility to run ce0l without physics is useless (ce0l is dedicated to Earth physics).
Location:
LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90

    r2331 r2336  
    11PROGRAM ce0l
    2   !
    3   ! Purpose: Calls etat0, creates initial states and limit_netcdf
    4   !
    5   ! interbar=.T. for barycentric interpolation inter_barxy
    6   ! extrap  =.T. for data extrapolation, like for the SSTs when file does not
    7   !                  contain ocean points only.
    8   ! oldice  =.T. for old-style ice, obtained using grille_m (grid_atob).
    9   ! masque is created in etat0, passed to limit to ensure consistancy.
    10 
    11   USE control_mod, only: DAY_STEP, DAYREF, NSPLIT_PHYS
    12   USE etat0dyn, only: etat0dyn_netcdf
    13   USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
     2!
     3!-------------------------------------------------------------------------------
     4! Purpose: Initial states and boundary conditions files creation:
     5!     * start.nc    for dynamics    (using etat0dyn     routine)
     6!     * startphy.nc for physics     (using etat0phys    routine)
     7!     * limit.nc    for forced runs (using limit_netcdf routine)
     8!-------------------------------------------------------------------------------
     9! Notes:
     10!     * extrap=.T. (default) for data extrapolation, like for the SSTs when file
     11!                   does contain ocean points only.
     12!     * "masque" can be:
     13!       - read from file "o2a.nc"          (for coupled runs).
     14!       - created in etat0phys or etat0dyn (for forced  runs).
     15!     It is then passed to limit_netcdf to ensure consistancy.
     16!-------------------------------------------------------------------------------
    1417  USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo
    15 
    16   USE etat0phys, only: etat0phys_netcdf
    17   USE dimphy, only: KLON
    18   USE infotrac, only: TYPE_TRAC, infotrac_init
     18  USE control_mod,    ONLY: day_step, dayref, nsplit_phys
     19  USE etat0dyn,       ONLY: etat0dyn_netcdf
     20  USE etat0phys,      ONLY: etat0phys_netcdf
     21  USE limit,          ONLY: limit_netcdf
     22  USE netcdf,         ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
     23  USE infotrac,       ONLY: type_trac, infotrac_init
     24  USE dimphy,         ONLY: klon
    1925  USE test_disvert_m, ONLY: test_disvert
     26  USE filtreg_mod,    ONLY: inifilr
     27#ifdef inca
     28  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
     29#endif
     30#ifdef CPP_PARA
     31  USE mod_const_mpi,  ONLY: init_const_mpi
     32  USE parallel_lmdz,  ONLY: init_parallel, mpi_rank, omp_rank, mpi_size
     33  USE bands,          ONLY: read_distrib, distrib_phys
     34  USE mod_hallo,      ONLY: init_mod_hallo
     35  USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
     36#endif
    2037
    2138  IMPLICIT NONE
    2239
    23   ! Local variables:
     40!-------------------------------------------------------------------------------
     41! Local variables:
    2442  include "dimensions.h"
    2543  include "paramet.h"
    26   include "comgeom.h"
     44  include "comgeom2.h"
    2745  include "comconst.h"
    2846  include "comvert.h"
     
    3048  include "temps.h"
    3149  include "logic.h"
    32   REAL               :: masque(iip1, jjp1)             !--- CONTINENTAL MASK
    33   REAL               :: phis  (iip1, jjp1)             !--- GROUND GEOPOTENTIAL
     50  REAL               :: masque(iip1,jjp1)             !--- CONTINENTAL MASK
     51  REAL               :: phis  (iip1,jjp1)             !--- GROUND GEOPOTENTIAL
    3452  CHARACTER(LEN=256) :: modname, fmt, calnd           !--- CALENDAR TYPE
    3553  LOGICAL            :: use_filtre_fft
    36   LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.
    37 
    38   !--- Local variables for ocean mask reading:
     54  LOGICAL, PARAMETER :: extrap=.FALSE.
     55
     56!--- Local variables for ocean mask reading:
    3957  INTEGER            :: nid_o2a, iml_omask, jml_omask, j
    4058  INTEGER            :: fid, iret, llm_tmp, ttm_tmp, itaul(1)
    41   REAL, ALLOCATABLE  :: lon_omask(:, :), dlon_omask(:), ocemask(:, :)
    42   REAL, ALLOCATABLE  :: lat_omask(:, :), dlat_omask(:), ocetmp (:, :)
     59  REAL, ALLOCATABLE  :: lon_omask(:,:), dlon_omask(:), ocemask(:,:)
     60  REAL, ALLOCATABLE  :: lat_omask(:,:), dlat_omask(:), ocetmp (:,:)
    4361  REAL               :: date, lev(1)
    44 
    45   !----------------------------------------------------------------------
     62!-------------------------------------------------------------------------------
    4663  modname="ce0l"
    4764
    48   !--- Constants
     65!--- Constants
    4966  pi     = 4. * ATAN(1.)
    5067  rad    = 6371229.
     
    5976
    6077  CALL conf_gcm( 99, .TRUE. )
    61 
    6278  dtvr = daysec/REAL(day_step)
    63   WRITE(lunout, *)'dtvr', dtvr
    64 
     79  WRITE(lunout,*)'dtvr',dtvr
    6580  CALL iniconst()
    6681  CALL inigeom()
    6782
     83!--- Calendar choice
    6884#ifdef CPP_IOIPSL
    6985  calnd='gregorian'
    7086  SELECT CASE(calend)
    71   CASE('earth_360d')
    72      CALL ioconf_calendar('360d')
    73      calnd='with 360 days/year'
    74   CASE('earth_365d')
    75      CALL ioconf_calendar('noleap')
    76      calnd='with no leap year'
    77   CASE('earth_366d')
    78      CALL ioconf_calendar('366d')
    79      calnd='with leap years only'
    80   CASE('gregorian')
    81      CALL ioconf_calendar('gregorian')
    82   CASE('standard')
    83      CALL ioconf_calendar('gregorian')
    84   CASE('julian')
    85      CALL ioconf_calendar('julian')
    86      calnd='julian'
    87   CASE('proleptic_gregorian')
    88      CALL ioconf_calendar('gregorian')
    89      !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
    90   CASE DEFAULT
    91      CALL abort_gcm('ce0l', 'Bad choice for calendar', 1)
     87    CASE('earth_360d');CALL ioconf_calendar('360d');   calnd='with 360 days/year'
     88    CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='with no leap year'
     89    CASE('earth_366d');CALL ioconf_calendar('366d');   calnd='with leap years only'
     90    CASE('gregorian'); CALL ioconf_calendar('gregorian')
     91    CASE('standard');  CALL ioconf_calendar('gregorian')
     92    CASE('julian');    CALL ioconf_calendar('julian'); calnd='julian'
     93    CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
     94  !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
     95    CASE DEFAULT
     96      CALL abort_gcm('ce0l','Bad choice for calendar',1)
    9297  END SELECT
    93   WRITE(lunout, *)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
    94 #endif
    95 
    96   use_filtre_fft=.FALSE.
    97   CALL getin('use_filtre_fft', use_filtre_fft)
    98   IF(use_filtre_fft) THEN
    99      WRITE(lunout, *)"FFT filter not available for sequential dynamics."
    100      WRITE(lunout, *)"Your setting of variable use_filtre_fft is not used."
    101   ENDIF
    102 
    103   !--- LAND MASK. TWO CASES:
    104   !   1) read from ocean model    file "o2a.nc"    (coupled runs)
    105   !   2) computed from topography file "Relief.nc" (masque(:, :)=-99999.)
    106   ! Coupled simulations (case 1) use the ocean model mask to compute the
    107   ! weights to ensure ocean fractions are the same for atmosphere and ocean.
    108 
    109   IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
    110      WRITE(lunout, *)'BEWARE !! No ocean mask "o2a.nc" file found'
    111      WRITE(lunout, *)'Forced run.'
    112      masque(:, :)=-99999.
    113   ELSE
    114      iret=NF90_CLOSE(nid_o2a)
    115      WRITE(lunout, *)'BEWARE !! Ocean mask "o2a.nc" file found'
    116      WRITE(lunout, *)'Coupled run.'
    117      CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
    118      IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
    119         WRITE(lunout, *)'Mismatching dimensions for ocean mask'
    120         WRITE(lunout, *)'iim  = ', iim , ' iml_omask = ', iml_omask
    121         WRITE(lunout, *)'jjp1 = ', jjp1, ' jml_omask = ', jml_omask
    122         CALL abort_gcm(modname, '', 1)
    123      END IF
    124      ALLOCATE(ocemask(iim, jjp1), lon_omask(iim, jjp1), dlon_omask(iim ))
    125      ALLOCATE(ocetmp (iim, jjp1), lat_omask(iim, jjp1), dlat_omask(jjp1))
    126      CALL flinopen("o2a.nc", .FALSE., iim, jjp1, llm_tmp, lon_omask, &
    127           lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
    128      CALL flinget(fid, "OceMask",    iim, jjp1, llm_tmp, ttm_tmp, 1, 1, ocetmp)
    129      CALL flinclo(fid)
    130      dlon_omask(1:iim ) = lon_omask(1:iim, 1)
    131      dlat_omask(1:jjp1) = lat_omask(1, 1:jjp1)
    132      ocemask = ocetmp
    133      IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
    134         DO j=1, jjp1
    135            ocemask(:, j) = ocetmp(:, jjp1-j+1)
    136         END DO
    137      END IF
    138      DEALLOCATE(ocetmp, lon_omask, lat_omask, dlon_omask, dlat_omask)
    139      IF(prt_level>=1) THEN
    140         WRITE(fmt, "(i4, 'i1)')")iim
    141         fmt='('//ADJUSTL(fmt)
    142         WRITE(lunout, *)'OCEAN MASK :'
    143         WRITE(lunout, fmt) NINT(ocemask)
    144      END IF
    145      masque(1:iim, :)=1.-ocemask(:, :)
    146      masque(iip1 , :)=masque(1, :)
    147      DEALLOCATE(ocemask)
    148   END IF
    149   phis(:, :)=-99999.
    150 
    151   CALL Init_Phys_lmdz(iim, jjp1, llm, 1, (/(jjm-1)*iim+2/))
    152   WRITE(lunout, *)'---> klon=', klon
    153 
    154   call infotrac_init
    155   CALL iniphysiq(iim, jjm, llm, daysec, dayref, dtphys / nsplit_phys, rlatu, &
    156        rlonv, aire, cu, cv, rad, g, r, cpp, iflag_phys)
    157 
    158   IF(pressure_exner) CALL test_disvert
    159 
     98  WRITE(lunout,*)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
     99#endif
     100
     101#ifdef CPP_PARA
     102!--- Physical grid + parallel initializations
     103  CALL init_const_mpi()
     104  CALL init_parallel()
     105  CALL read_distrib()
     106  CALL init_mod_hallo()
     107  CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
     108  CALL init_interface_dyn_phys()
     109#else
     110  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     111#endif
     112  WRITE(lunout,*)'---> klon=',klon
     113
     114!--- Tracers initializations
    160115  IF (type_trac == 'inca') THEN
    161116#ifdef INCA
    162      CALL init_const_lmdz(nbtr, anneeref, dayref, iphysiq, day_step, nday)
    163      CALL init_inca_para(iim, jjm+1, klon, 1, klon_mpi_para_nb, 0)
    164      WRITE(lunout, *)'nbtr =' , nbtr
    165 #endif
    166   END IF
     117    CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,&
     118                         nbsrf,is_oce,is_sic,is_ter,is_lic,calend)
     119    CALL init_inca_para(iim,jjp1,llm,klon_glo,mpi_size,distrib_phys,&
     120                        COMM_LMDZ)
     121    WRITE(lunout,*)'nbtr =' , nbtr
     122#endif
     123  END IF
     124  CALL infotrac_init()
     125
     126  CALL inifilr()
     127  CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, &
     128                 rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
     129  IF(pressure_exner) CALL test_disvert
     130
     131#ifdef CPP_PARA
     132  IF (mpi_rank==0.AND.omp_rank==0) THEN
     133#endif
     134  use_filtre_fft=.FALSE.
     135  CALL getin('use_filtre_fft',use_filtre_fft)
     136  IF(use_filtre_fft) THEN
     137     WRITE(lunout,*)"FFT filter not available for sequential dynamics."
     138     WRITE(lunout,*)"Your setting of variable use_filtre_fft is not used."
     139  ENDIF
     140
     141!--- LAND MASK. TWO CASES:
     142!   1) read from ocean model    file "o2a.nc"    (coupled runs)
     143!   2) computed from topography file "Relief.nc" (masque(:,:)=-99999.)
     144! Coupled simulations (case 1) use the ocean model mask to compute the
     145! weights to ensure ocean fractions are the same for atmosphere and ocean.
     146!*******************************************************************************
     147  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
     148    WRITE(lunout,*)'BEWARE !! No ocean mask "o2a.nc" file found'
     149    WRITE(lunout,*)'Forced run.'
     150    masque(:,:)=-99999.
     151  ELSE
     152    iret=NF90_CLOSE(nid_o2a)
     153    WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found'
     154    WRITE(lunout,*)'Coupled run.'
     155    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
     156    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
     157      WRITE(lunout,*)'Mismatching dimensions for ocean mask'
     158      WRITE(lunout,*)'iim  = ',iim ,' iml_omask = ',iml_omask
     159      WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
     160      CALL abort_gcm(modname,'',1)
     161    END IF
     162    ALLOCATE(ocemask(iim,jjp1),lon_omask(iim,jjp1),dlon_omask(iim ))
     163    ALLOCATE(ocetmp (iim,jjp1),lat_omask(iim,jjp1),dlat_omask(jjp1))
     164    CALL flinopen("o2a.nc", .FALSE.,iim,jjp1,llm_tmp,lon_omask,lat_omask,lev,  &
     165     &                                       ttm_tmp,itaul,date,dt,fid)
     166    CALL flinget(fid, "OceMask",    iim,jjp1,llm_tmp,ttm_tmp,1,1,ocetmp)
     167    CALL flinclo(fid)
     168    dlon_omask(1:iim ) = lon_omask(1:iim,1)
     169    dlat_omask(1:jjp1) = lat_omask(1,1:jjp1)
     170    ocemask = ocetmp
     171    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
     172      DO j=1,jjp1; ocemask(:,j) = ocetmp(:,jjp1-j+1); END DO
     173    END IF
     174    DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask)
     175    IF(prt_level>=1) THEN
     176      WRITE(fmt,"(i4,'i1)')")iim ; fmt='('//ADJUSTL(fmt)
     177      WRITE(lunout,*)'OCEAN MASK :'
     178      WRITE(lunout,fmt) NINT(ocemask)
     179    END IF
     180    masque(1:iim,:)=1.-ocemask(:,:)
     181    masque(iip1 ,:)=masque(1,:)
     182    DEALLOCATE(ocemask)
     183  END IF
     184  phis(:,:)=-99999.
     185
    167186  IF(ok_etat0) THEN
    168      WRITE(lunout, '(//)')
    169      WRITE(lunout, *) '  ************************  '
    170      WRITE(lunout, *) '  ***  etat0phy_netcdf ***  '
    171      WRITE(lunout, *) '  ************************  '
    172      WRITE(lunout, '(//)')
    173      WRITE(lunout, *) ' interbar = ', interbar
    174      CALL etat0phys_netcdf(interbar, masque, phis)
    175   END IF
    176 
    177   IF(ok_etat0) THEN
    178      WRITE(lunout, '(//)')
    179      WRITE(lunout, *) '  ************************  '
    180      WRITE(lunout, *) '  ***  etat0dyn_netcdf ***  '
    181      WRITE(lunout, *) '  ************************  '
    182      WRITE(lunout, '(//)')
    183      WRITE(lunout, *) ' interbar = ', interbar
    184      CALL etat0dyn_netcdf(interbar, masque, phis)
     187    WRITE(lunout,'(//)')
     188    WRITE(lunout,*) '  ************************  '
     189    WRITE(lunout,*) '  ***  etat0phy_netcdf ***  '
     190    WRITE(lunout,*) '  ************************  '
     191    CALL etat0phys_netcdf(masque,phis)
     192    WRITE(lunout,'(//)')
     193    WRITE(lunout,*) '  ************************  '
     194    WRITE(lunout,*) '  ***  etat0dyn_netcdf ***  '
     195    WRITE(lunout,*) '  ************************  '
     196    CALL etat0dyn_netcdf(masque,phis)
    185197  END IF
    186198
    187199  IF(ok_limit) THEN
    188      WRITE(lunout, '(//)')
    189      WRITE(lunout, *) '  *********************  '
    190      WRITE(lunout, *) '  ***  Limit_netcdf ***  '
    191      WRITE(lunout, *) '  *********************  '
    192      WRITE(lunout, '(//)')
    193      CALL limit_netcdf(interbar, extrap, oldice, masque)
    194   END IF
    195 
    196   WRITE(lunout, '(//)')
    197   WRITE(lunout, *) '  ***************************  '
    198   WRITE(lunout, *) '  ***  grilles_gcm_netcdf ***  '
    199   WRITE(lunout, *) '  ***************************  '
    200   WRITE(lunout, '(//)')
    201   CALL grilles_gcm_netcdf_sub(masque, phis)
     200    WRITE(lunout,'(//)')
     201    WRITE(lunout,*) '  *********************  '
     202    WRITE(lunout,*) '  ***  Limit_netcdf ***  '
     203    WRITE(lunout,*) '  *********************  '
     204    WRITE(lunout,'(//)')
     205    CALL limit_netcdf(masque,phis,extrap)
     206  END IF
     207
     208  WRITE(lunout,'(//)')
     209  WRITE(lunout,*) '  ***************************  '
     210  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
     211  WRITE(lunout,*) '  ***************************  '
     212  WRITE(lunout,'(//)')
     213  CALL grilles_gcm_netcdf_sub(masque,phis)
     214
     215#ifdef CPP_PARA
     216  END IF
     217#endif
    202218
    203219END PROGRAM ce0l
     220!
     221!-------------------------------------------------------------------------------
  • LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90

    r2334 r2336  
    1212!  routine (to be called after restget):
    1313!    CALL startget_dyn3d(varname, lon_in,  lat_in, pls, workvar,&
    14 !                          champ, lon_in2, lat_in2, ibar)
     14!                          champ, lon_in2, lat_in2)
    1515!
    1616!    *  Variables should have the following names in the NetCDF files:
     
    3636  USE ioipsl,         ONLY: flininfo, flinopen, flinget, flinclo, histclo
    3737  USE assert_eq_m,    ONLY: assert_eq
    38 #ifdef CPP_PHYS
    3938  USE indice_sol_mod, ONLY: epsfra
    40 #endif
    4139  IMPLICIT NONE
    4240
     
    5452  include "serre.h"
    5553  REAL, SAVE :: deg2rad
    56 #ifndef CPP_PHYS
    57   REAL, SAVE :: epsfra= 1.E-5
    58 #endif
    5954  INTEGER,            SAVE      :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
    6055  REAL, ALLOCATABLE,  SAVE      :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
    6156  CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
    62   CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc'
    6357
    6458CONTAINS
     
    6660!-------------------------------------------------------------------------------
    6761!
    68 SUBROUTINE etat0dyn_netcdf(ib, masque, phis)
     62SUBROUTINE etat0dyn_netcdf(masque, phis)
    6963!
    7064!-------------------------------------------------------------------------------
     
    7367! Notes:  1) This routine is designed to work for Earth
    7468!         2) If masque(:,:)/=-99999., masque and phis are already known.
    75 !         Otherwise: read it from ocean model file (coupled run) or compute it.
     69!         Otherwise: compute it.
    7670!-------------------------------------------------------------------------------
    7771  USE control_mod
    78 #ifdef CPP_PHYS
    7972  USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
    8073  USE regr_pr_o3_m,   ONLY: regr_pr_o3
    8174  USE press_coefoz_m, ONLY: press_coefoz
    82 #endif
    8375  USE exner_hyb_m,    ONLY: exner_hyb
    8476  USE exner_milieu_m, ONLY: exner_milieu
    85   USE infotrac, only: NQTOT, TNAME
     77  USE infotrac,       ONLY: nqtot, tname
    8678  USE filtreg_mod
    8779  IMPLICIT NONE
    8880!-------------------------------------------------------------------------------
    8981! Arguments:
    90   LOGICAL, INTENT(IN)    :: ib                  !--- Barycentric interpolation
    9182  REAL,    INTENT(INOUT) :: masque(iip1,jjp1)   !--- Land-ocean mask
    9283  REAL,    INTENT(INOUT) :: phis  (iip1,jjp1)   !--- Ground geopotential
     
    111102  deg2rad = pi/180.0
    112103
    113 ! Initializations for tracers and filter
    114 !*******************************************************************************
    115   CALL inifilr()
    116 
    117104! Compute ground geopotential and possibly the mask.
    118105!*******************************************************************************
    119106  masque_tmp(:,:)=masque(:,:)
    120   CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp)
    121107  WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
    122108  IF(ALL(masque==-99999.)) THEN                         !--- KEEP NEW MASK
     
    132118! Compute psol AND tsol, knowing phis.
    133119!*******************************************************************************
    134   CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, ib, phis, psol)
     120  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol)
    135121
    136122! Mid-levels pressure computation
     
    147133!*******************************************************************************
    148134  uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
    149   CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv,ib)
     135  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv)
    150136  CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
    151  &                           rlonu,rlatu(:jjm),ib)
    152   CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv,ib)
     137 &                           rlonu,rlatu(:jjm))
     138  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv)
    153139  tpot(:,:,:)=t3d(:,:,:)
    154   CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv,ib)
     140  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv)
    155141
    156142  WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
     
    165151!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
    166152  qd (:,:,:) = 0.0
    167   CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv,ib)
     153  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv)
    168154  ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
    169155  CALL flinclo(fid_dyn)
    170156
    171 #ifdef CPP_PHYS
    172157! Parameterization of ozone chemistry:
    173158!*******************************************************************************
     
    180165    q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
    181166  END IF
    182 
    183 #endif
    184167  q3d(iip1,:,:,:)=q3d(1,:,:,:)
    185 
    186 ! Intermediate computation
    187 !*******************************************************************************
    188   CALL massdair(p3d,masse)
    189   WRITE(lunout,*)' ALPHAX ',alphax
    190   DO l=1,llm
    191     xppn(:)=aire(1:iim,1   )*masse(1:iim,1   ,l)
    192     xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
    193     xpn=SUM(xppn)/apoln
    194     xps=SUM(xpps)/apols
    195     masse(:,1   ,l)=xpn
    196     masse(:,jjp1,l)=xps
    197   END DO
    198168
    199169! Writing
     
    215185  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
    216186  WRITE(lunout,*)'sortie geopot'
    217   CALL caldyn0( itau, uvent, vvent, tpot, psol, pk, phis,               &
    218                 phi,  w, time+iday-dayref)
     187  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
     188                phi,  w, pbaru, pbarv, time+iday-dayref)
    219189  WRITE(lunout,*)'sortie caldyn0'     
     190#ifdef CPP_PARA
     191  CALL dynredem0_loc( "start.nc", dayref, phis)
     192#else
    220193  CALL dynredem0( "start.nc", dayref, phis)
     194#endif
    221195  WRITE(lunout,*)'sortie dynredem0'
     196#ifdef CPP_PARA
     197  CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
     198#else
    222199  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
     200#endif
    223201  WRITE(lunout,*)'sortie dynredem1'
    224202  CALL histclo()
     
    232210!
    233211SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
    234                         champ, lon_in2, lat_in2, ibar)
     212                        champ, lon_in2, lat_in2)
    235213!-------------------------------------------------------------------------------
    236214  IMPLICIT NONE
     
    252230  REAL,             INTENT(IN)    :: lon_in2(:)       ! dim (iml)
    253231  REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
    254   LOGICAL,          INTENT(IN)    :: ibar
    255232!-------------------------------------------------------------------------------
    256233! Local variables:
     
    287264!--- INTERPOLATE 3D FIELD IF NEEDED
    288265  IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
    289                                                   lat_in2,pls,champ,ibar)
     266                                                  lat_in2,pls,champ)
    290267
    291268!--- COMPUTE THE REQUIRED FILED
     
    317294!-------------------------------------------------------------------------------
    318295!
    319 SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
    320 !
    321 !-------------------------------------------------------------------------------
    322   USE conf_dat_m, ONLY: conf_dat2d
    323   IMPLICIT NONE
    324 !===============================================================================
    325 ! Purpose:  Compute "phis" just like it would be in start_init_orog.
    326 !===============================================================================
    327 ! Arguments:
    328   REAL, INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
    329   REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
    330 !-------------------------------------------------------------------------------
    331 ! Local variables:
    332   CHARACTER(LEN=256) :: modname="start_init_orog0"
    333   CHARACTER(LEN=256) :: title="RELIEF"
    334   INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
    335   REAL               :: lev(1), date, dt
    336   REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
    337   REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
    338 !-------------------------------------------------------------------------------
    339   iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
    340   jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
    341   IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
    342   IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
    343   pi=2.0*ASIN(1.0); deg2rad=pi/180.0
    344   IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
    345 
    346 !--- HIGH RESOLUTION OROGRAPHY
    347   CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
    348 
    349   ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
    350   CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
    351                 lev, ttm_tmp, itau, date, dt, fid)
    352   ALLOCATE(relief_hi(iml_rel,jml_rel))
    353   CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
    354   CALL flinclo(fid)
    355 
    356 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
    357   ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
    358   lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
    359   lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
    360 
    361 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
    362   ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
    363   CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
    364   DEALLOCATE(lon_ini,lat_ini)
    365 
    366 !--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
    367   WRITE(lunout,*)
    368   WRITE(lunout,*)'*** Compute surface geopotential ***'
    369 
    370 !--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
    371   CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
    372   phis = phis * 9.81
    373   phis(iml,:) = phis(1,:)
    374   DEALLOCATE(relief_hi,lon_rad,lat_rad)
    375 
    376 END SUBROUTINE start_init_orog0
    377 !
    378 !-------------------------------------------------------------------------------
    379 
    380 
    381 !-------------------------------------------------------------------------------
    382 !
    383 SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
    384 !
    385 !===============================================================================
    386 ! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
    387 !          without any call to physics subroutines.
    388 !===============================================================================
    389   IMPLICIT NONE
    390 !-------------------------------------------------------------------------------
    391 ! Arguments:
    392   REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
    393   REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
    394   REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
    395   REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
    396   REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
    397 !-------------------------------------------------------------------------------
    398 ! Local variables:
    399   CHARACTER(LEN=256) :: modname="grid_noro0"
    400   REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
    401   REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
    402   REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
    403   REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
    404   REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
    405   REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
    406   REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
    407   LOGICAL :: masque_lu
    408   INTEGER :: i, ii, imdp, imar, iext
    409   INTEGER :: j, jj, jmdp, jmar, nn
    410   REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
    411   REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
    412 !-------------------------------------------------------------------------------
    413   imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
    414   jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
    415   imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
    416   jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
    417   IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
    418   IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
    419   iext=imdp/10
    420   xpi = ACOS(-1.)
    421   rad = 6371229.
    422 
    423 !--- ARE WE USING A READ MASK ?
    424   masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
    425   WRITE(lunout,*)'Masque lu: ',masque_lu
    426 
    427 !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
    428   ALLOCATE(xusn(imdp+2*iext))
    429   xusn(1     +iext:imdp  +iext)=xd(:)
    430   xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
    431   xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
    432 
    433   ALLOCATE(yusn(jmdp+2))
    434   yusn(1       )=yd(1)   +(yd(1)   -yd(2))
    435   yusn(2:jmdp+1)=yd(:)
    436   yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
    437 
    438   ALLOCATE(zusn(imdp+2*iext,jmdp+2))
    439   zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
    440   zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
    441   zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
    442   zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
    443   zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
    444   zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
    445   zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
    446 
    447 !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
    448   ALLOCATE(a(imar+1),b(imar+1))
    449   b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
    450   b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
    451   a(1)=x(1)-(x(2)-x(1))/2.0
    452   a(2:imar+1)= b(1:imar)
    453 
    454   ALLOCATE(c(jmar),d(jmar))
    455   d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
    456   d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
    457   c(1)=y(1)-(y(2)-y(1))/2.0
    458   c(2:jmar)=d(1:jmar-1)
    459 
    460 !--- INITIALIZATIONS:
    461   ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
    462   ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
    463 
    464 !--- SUMMATION OVER GRIDPOINT AREA
    465   zleny=xpi/REAL(jmdp)*rad
    466   xincr=xpi/REAL(jmdp)/2.
    467   ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
    468   ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
    469   DO ii = 1, imar+1
    470     DO jj = 1, jmar
    471       DO j = 2,jmdp+1
    472         zlenx  =zleny  *COS(yusn(j))
    473         zbordnor=(xincr+c(jj)-yusn(j))*rad
    474         zbordsud=(xincr-d(jj)+yusn(j))*rad
    475         weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
    476         IF(weighy/=0) THEN
    477           DO i = 2, imdp+2*iext-1
    478             zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
    479             zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
    480             weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
    481             IF(weighx/=0)THEN
    482               num_tot(ii,jj)=num_tot(ii,jj)+1.0
    483               IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
    484               weight(ii,jj)=weight(ii,jj)+weighx*weighy
    485               zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
    486             END IF
    487           END DO
    488         END IF
    489       END DO
    490     END DO
    491   END DO
    492 
    493 !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
    494   IF(.NOT.masque_lu) THEN
    495     WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
    496   END IF
    497   nn=COUNT(weight(:,1:jmar-1)==0.0)
    498   IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
    499   WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
    500 
    501 !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
    502   ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
    503   WHERE(mask>=0.1) mask_tmp = 1.
    504   WHERE(weight(:,:)/=0.0)
    505     zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
    506     zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
    507   END WHERE
    508   WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
    509 
    510 !--- Values at poles
    511   zphi(imar+1,:)=zphi(1,:)
    512 
    513   zweinor=SUM(weight(1:imar,   1),DIM=1)
    514   zweisud=SUM(weight(1:imar,jmar),DIM=1)
    515   zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
    516   zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
    517   zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
    518 
    519 END SUBROUTINE grid_noro0
    520 !
    521 !-------------------------------------------------------------------------------
    522 
    523 
    524 !-------------------------------------------------------------------------------
    525 !
    526 SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol)
     296SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol)
    527297!
    528298!-------------------------------------------------------------------------------
     
    534304  REAL,    INTENT(IN)  :: lon_in (:),  lat_in (:)    ! dim (iml) (jml)
    535305  REAL,    INTENT(IN)  :: lon_in2(:),  lat_in2(:)    ! dim (iml) (jml2)
    536   LOGICAL, INTENT(IN)  :: ibar
    537306  REAL,    INTENT(IN)  :: zs  (:,:)                  ! dim (iml,jml)
    538307  REAL,    INTENT(OUT) :: psol(:,:)                  ! dim (iml,jml)
     
    606375    ALLOCATE(field(iml,jml))
    607376    CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
    608     CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)
    609     CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,       &
    610                               lon_in, lat_in, lon_in2, lat_in2, field)
     377    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
     378    CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana,              &
     379                                        lon_in, lat_in, lon_in2, lat_in2, field)
    611380  ELSE IF(SIZE(field)/=SIZE(z)) THEN
    612381    msg='The '//TRIM(msg)//' field we have does not have the right size'
     
    625394!-------------------------------------------------------------------------------
    626395!
    627 SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d,ibar)
     396SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d)
    628397!
    629398!-------------------------------------------------------------------------------
     
    639408  REAL,    INTENT(IN)  :: pls_in(:,:,:)           ! dim (iml,jml,lml)
    640409  REAL,    INTENT(OUT) :: var3d (:,:,:)           ! dim (iml,jml,lml)
    641   LOGICAL, INTENT(IN)  :: ibar
    642410!-------------------------------------------------------------------------------
    643411! Local variables:
     
    667435  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
    668436  CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini,                           &
    669                        lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
     437                       lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.)
    670438  DEALLOCATE(lon_ini, lat_ini)
    671439
     
    673441  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
    674442  DO il = 1,llm_dyn
    675     CALL interp_startvar(var,ibar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),     &
     443    CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),          &
    676444                          lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
    677445  END DO
     
    708476!-------------------------------------------------------------------------------
    709477!
    710 SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
     478SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
    711479!
    712480!-------------------------------------------------------------------------------
    713481  USE inter_barxy_m, ONLY: inter_barxy
    714   USE grid_atob_m,   ONLY: grille_m
    715482  IMPLICIT NONE
    716483!-------------------------------------------------------------------------------
    717484! Arguments:
    718485  CHARACTER(LEN=*), INTENT(IN)  :: nam
    719   LOGICAL,          INTENT(IN)  :: ibar, ibeg
     486  LOGICAL,          INTENT(IN)  :: ibeg
    720487  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
    721488  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
     
    735502  j2=SIZE(lat2)
    736503  ALLOCATE(vtmp(i1-1,j1))
    737   IF(ibar) THEN
    738     IF(ibeg.AND.prt_level>1) THEN
    739       WRITE(lunout,*)"---------------------------------------------------------"
    740       WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
    741       WRITE(lunout,*)"---------------------------------------------------------"
    742     END IF
    743     CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
    744   ELSE
    745     CALL grille_m   (lon, lat,        vari, lon1,        lat1, vtmp)
    746   END IF
     504  IF(ibeg.AND.prt_level>1) THEN
     505    WRITE(lunout,*)"---------------------------------------------------------"
     506    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
     507    WRITE(lunout,*)"---------------------------------------------------------"
     508  END IF
     509  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
    747510  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
    748511
  • LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90

    r2320 r2336  
    3636  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
    3737  USE assert_eq_m,        ONLY: assert_eq
    38   USE dimphy
     38  USE dimphy,             ONLY: klon
     39  USE conf_dat_m,         ONLY: conf_dat2d
    3940  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
    4041    rlon, solsw, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
     
    5152  include "paramet.h"
    5253  include "comgeom2.h"
    53   include "comvert.h"
    5454  include "comconst.h"
    5555  include "dimsoil.h"
    5656  include "temps.h"
    57   include "comdissnew.h"
    58   include "serre.h"
    5957  include "clesphys.h"
    6058  REAL, SAVE :: deg2rad
    6159  REAL, SAVE, ALLOCATABLE :: tsol(:)
    62 !  REAL, SAVE, ALLOCATABLE :: rugo(:,:)  ! ??? COMPUTED BUT NOT USED ???
    6360  INTEGER,            SAVE      :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
    6461  REAL, ALLOCATABLE,  SAVE      :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
    65   CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc"
    66   CHARACTER(LEN=256), PARAMETER :: title="RELIEF"
     62  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
     63  CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc",  psrfvar="SP"
     64  CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW",       tsrfvar="ST"
    6765
    6866
     
    7270!-------------------------------------------------------------------------------
    7371!
    74 SUBROUTINE etat0phys_netcdf(ib, masque, phis)
     72SUBROUTINE etat0phys_netcdf(masque, phis)
    7573!
    7674!-------------------------------------------------------------------------------
    7775! Purpose: Creates initial states
    7876!-------------------------------------------------------------------------------
    79 ! Note: This routine is designed to work for Earth
     77! Notes:  1) This routine is designed to work for Earth
     78!         2) If masque(:,:)/=-99999., masque and phis are already known.
     79!         Otherwise: compute it.
    8080!-------------------------------------------------------------------------------
    8181  USE control_mod
     
    8484  USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
    8585  USE indice_sol_mod
    86   USE conf_phys_m,    ONLY: conf_phys
    87   USE exner_hyb_m,    ONLY: exner_hyb
    88   USE exner_milieu_m, ONLY: exner_milieu
    89   USE test_disvert_m, ONLY: test_disvert
    90   USE grid_atob_m,    ONLY: grille_m
     86  USE conf_phys_m, ONLY: conf_phys
     87  USE init_ssrf_m, ONLY: start_init_subsurf
    9188  IMPLICIT NONE
    9289!-------------------------------------------------------------------------------
    9390! Arguments:
    94   LOGICAL, INTENT(IN)    :: ib          !--- Barycentric interpolation
    9591  REAL,    INTENT(INOUT) :: masque(:,:) !--- Land mask           dim(iip1,jjp1)
    9692  REAL,    INTENT(INOUT) :: phis  (:,:) !--- Ground geopotential dim(iip1,jjp1)
     
    9995  CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
    10096  INTEGER            :: i, j, l, ji, iml, jml
    101   REAL               :: phystep
     97  LOGICAL            :: read_mask
     98  REAL               :: phystep, dummy
    10299  REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp
    103100  REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
    104101  REAL, DIMENSION(klon,nbsrf)         :: qsolsrf, snsrf
    105102  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
    106 
    107 !--- Local variables for sea-ice reading:
    108   LOGICAL           :: read_mask
    109   INTEGER           :: iml_lic, jml_lic, isst(klon-2)
    110   INTEGER           :: fid, llm_tmp, ttm_tmp, itaul(1)
    111   REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:)
    112   REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:)
    113   REAL              :: date, lev(1), dummy
    114   REAL              :: flic_tmp(SIZE(masque,1),SIZE(masque,2))
    115103
    116104!--- Arguments for conf_phys
     
    132120  jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
    133121
    134 ! Grid construction and miscellanous initializations.
     122! Physics configuration
    135123!*******************************************************************************
    136124  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
     
    145133                   read_climoz,                                         &
    146134                   alp_offset)
    147 
    148135  CALL phys_state_var_init(read_climoz)
    149136
    150   co2_ppm0 = co2_ppm  !--- Initial atmospheric CO2 conc. from .def file
    151 
    152   rlat(1) = pi/2.
    153   DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);     END DO
    154   rlat(klon) = - pi/2.
    155   rlat(:)=rlat(:)*(180.0/pi)
    156 
    157   rlon(1) = 0.0
    158   DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO
    159   rlon(klon) = 0.0
    160   rlon(:)=rlon(:)*(180.0/pi)
     137!--- Initial atmospheric CO2 conc. from .def file
     138  co2_ppm0 = co2_ppm
    161139
    162140! Compute ground geopotential, sub-cells quantities and possibly the mask.
     
    174152    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
    175153  END IF
    176  CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
     154  CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
    177155
    178156! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
    179157!*******************************************************************************
    180   CALL start_init_phys(rlonv, rlatu, rlonu, rlatv, ib, phis)
     158  CALL start_init_phys(rlonu, rlatv, phis)
    181159
    182160! Some initializations.
     
    188166    CALL regr_lat_time_climoz(read_climoz)
    189167
    190 ! Sub-surfaces initialization
    191 !*******************************************************************************
    192 !--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
    193   CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
    194   ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
    195   ALLOCATE(dlat_lic(jml_lic),       dlon_lic(iml_lic))
    196   ALLOCATE( fraclic(iml_lic,jml_lic))
    197   CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp,  &
    198  &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
    199   CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
    200   CALL flinclo(fid)
    201   WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic
    202   IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees
    203   IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
    204   dlon_lic(:)=lon_lic(:,1)
    205   dlat_lic(:)=lat_lic(1,:)
    206   CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) )
    207   flic_tmp(iml,:)=flic_tmp(1,:)
    208 
    209 !--- To the physical grid
    210   pctsrf(:,:) = 0.
    211   CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic))
    212 
    213 !--- Adequation with soil/sea mask
    214   WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0.
    215   WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
    216   pctsrf(:,is_ter)=zmasq(:)
    217   DO ji=1,klon
    218     IF(zmasq(ji)>EPSFRA) THEN
    219       IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
    220         pctsrf(ji,is_lic)=zmasq(ji)
    221         pctsrf(ji,is_ter)=0.
    222       ELSE
    223         pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
    224         IF(pctsrf(ji,is_ter)<EPSFRA) THEN
    225           pctsrf(ji,is_ter)=0.
    226           pctsrf(ji,is_lic)=zmasq(ji)
    227         END IF
    228       END IF
    229     END IF
    230   END DO
    231 
    232 !--- Sub-surface ocean and sea ice (sea ice set to zero for start).
    233   pctsrf(:,is_oce)=(1.-zmasq(:))
    234   WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
    235   IF(read_mask) pctsrf(:,is_oce)=1-zmasq(:)
    236   isst=0
    237   WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
    238 
    239 !--- It is checked that the sub-surfaces sum is equal to 1.
    240   ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
    241   IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points'
     168! Sub-surfaces initialization.
     169!*******************************************************************************
     170  CALL start_init_subsurf(read_mask)
    242171
    243172! Write physical initial state
     
    321250!   described in LOTT & MILLER (1997) and LOTT(1999).
    322251!===============================================================================
    323   USE conf_dat_m,  ONLY: conf_dat2d
    324 !  USE grid_atob_m, ONLY: rugsoro
    325252  USE grid_noro_m, ONLY: grid_noro
    326253  IMPLICIT NONE
     
    331258!-------------------------------------------------------------------------------
    332259! Local variables:
    333   CHARACTER(LEN=256) :: modname="start_init_orog"
    334   CHARACTER(LEN=256) :: title="RELIEF"
     260  CHARACTER(LEN=256) :: modname
    335261  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
    336262  REAL               :: lev(1), date, dt
     
    340266  REAL, ALLOCATABLE  :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
    341267!-------------------------------------------------------------------------------
     268  modname="start_init_orog"
    342269  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
    343270  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
     
    350277                lev, ttm_tmp, itau, date, dt, fid)
    351278  ALLOCATE(relief_hi(iml_rel,jml_rel))
    352   CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
     279  CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)
    353280  CALL flinclo(fid)
    354281
     
    360287!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
    361288  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
    362   CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
     289  CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)
    363290  DEALLOCATE(lon_ini,lat_ini)
    364291
     
    378305  phis = phis * 9.81
    379306  phis(iml,:) = phis(1,:)
    380 
    381 !--- COMPUTE SURFACE ROUGHNESS
    382 !  WRITE(lunout,*)
    383 !  WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***'
    384 !  ALLOCATE(tmp_var(iml-1,jml))
    385 !  CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var)
    386 !  ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
    387 !  DEALLOCATE(tmp_var)
    388307  DEALLOCATE(relief_hi,lon_rad,lat_rad)
    389308
     
    405324!-------------------------------------------------------------------------------
    406325!
    407 SUBROUTINE start_init_phys(lon_in,lat_in,lon_in2,lat_in2,ibar,phis)
     326SUBROUTINE start_init_phys(lon_in,lat_in,phis)
    408327!
    409328!===============================================================================
     
    413332!-------------------------------------------------------------------------------
    414333! Arguments:
    415   REAL,    INTENT(IN) :: lon_in (:),  lat_in (:)     ! dim (iml) (jml)
    416   REAL,    INTENT(IN) :: lon_in2(:),  lat_in2(:)     ! dim (iml) (jml2)
    417   LOGICAL, INTENT(IN) :: ibar
     334  REAL,    INTENT(IN) :: lon_in(:),  lat_in(:)       ! dim (iml) (jml2)
    418335  REAL,    INTENT(IN) :: phis(:,:)                   ! dim (iml,jml)
    419336!-------------------------------------------------------------------------------
    420337! Local variables:
    421   CHARACTER(LEN=256) :: modname="start_init_phys", physfname="ECPHY.nc"
     338  CHARACTER(LEN=256) :: modname
    422339  REAL               :: date, dt
    423340  INTEGER            :: iml, jml, jml2, itau(1)
     
    426343  REAL, ALLOCATABLE  :: ts(:,:), qs(:,:)
    427344!-------------------------------------------------------------------------------
    428   iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml")
    429   jml=assert_eq(SIZE(lat_in),SIZE(phis,2),              TRIM(modname)//" jml")
    430   jml2=SIZE(lat_in2)
     345  modname="start_init_phys"
     346  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
     347  jml=SIZE(phis,2); jml2=SIZE(lat_in)
    431348
    432349  WRITE(lunout,*)'Opening the surface analysis'
    433   CALL flininfo(physfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
    434   WRITE(lunout,*) 'Values read: ',   iml_phys, jml_phys, llm_phys, ttm_phys
     350  CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
     351  WRITE(lunout,*) 'Values read: ',  iml_phys, jml_phys, llm_phys, ttm_phys
    435352
    436353  ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
    437354  ALLOCATE(levphys_ini(llm_phys))
    438   CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
     355  CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
    439356                lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
    440357
     
    445362
    446363  ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
    447   CALL get_var_phys('ST'  ,ts)                   !--- SURFACE TEMPERATURE
    448   CALL get_var_phys('CDSW',qs)                   !--- SOIL MOISTURE
     364  CALL get_var_phys(tsrfvar,ts)                   !--- SURFACE TEMPERATURE
     365  CALL get_var_phys(qsolvar,qs)                   !--- SOIL MOISTURE
    449366  CALL flinclo(fid_phys)
    450367  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
     
    463380!
    464381!-------------------------------------------------------------------------------
    465   USE conf_dat_m, ONLY: conf_dat2d
    466382  IMPLICIT NONE
    467383!-------------------------------------------------------------------------------
     
    474390!-------------------------------------------------------------------------------
    475391  SELECT CASE(title)
    476     CASE('SP');        tllm=0
    477     CASE('ST','CDSW'); tllm=llm_phys
     392    CASE(psrfvar);         tllm=0
     393    CASE(tsrfvar,qsolvar); tllm=llm_phys
    478394  END SELECT
    479395  IF(ALLOCATED(field)) RETURN
    480396  ALLOCATE(field(iml,jml)); field(:,:)=0.
    481397  CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana)
    482   CALL conf_dat2d(title, lon_ini, lat_ini,  lon_rad, lat_rad, var_ana, ibar)
    483   CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,         &
    484                             lon_in, lat_in, lon_in2, lat_in2, field)
     398  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
     399  CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana,               &
     400                                      lon_in,  lat_in, field)
    485401
    486402END SUBROUTINE get_var_phys
     
    495411!-------------------------------------------------------------------------------
    496412!
    497 SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
     413SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
    498414!
    499415!-------------------------------------------------------------------------------
    500416  USE inter_barxy_m, ONLY: inter_barxy
    501   USE grid_atob_m,   ONLY: grille_m
    502417  IMPLICIT NONE
    503418!-------------------------------------------------------------------------------
    504419! Arguments:
    505420  CHARACTER(LEN=*), INTENT(IN)  :: nam
    506   LOGICAL,          INTENT(IN)  :: ibar, ibeg
     421  LOGICAL,          INTENT(IN)  :: ibeg
    507422  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
    508423  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
    509   REAL,             INTENT(IN)  :: lon1(:), lat1(:) ! dim (i1) (j1)
    510424  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
    511425  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
    512426!-------------------------------------------------------------------------------
    513427! Local variables:
    514   CHARACTER(LEN=256) :: modname="interp_startvar"
     428  CHARACTER(LEN=256) :: modname
    515429  INTEGER            :: ii, jj, i1, j1, j2
    516430  REAL, ALLOCATABLE  :: vtmp(:,:)
    517431!-------------------------------------------------------------------------------
    518   ii=assert_eq(SIZE(lon),            SIZE(vari,1),TRIM(modname)//" ii")
    519   jj=assert_eq(SIZE(lat),            SIZE(vari,2),TRIM(modname)//" jj")
    520   i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
    521   j1=assert_eq(SIZE(lat1),           SIZE(varo,2),TRIM(modname)//" j1")
    522   j2=SIZE(lat2)
     432  modname="interp_startvar"
     433  ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
     434  jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
     435  i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
     436  j1=SIZE(varo,2); j2=SIZE(lat2)
    523437  ALLOCATE(vtmp(i1-1,j1))
    524   IF(ibar) THEN
    525     IF(ibeg.AND.prt_level>1) THEN
    526       WRITE(lunout,*)"--------------------------------------------------------"
    527       WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
    528       WRITE(lunout,*)"--------------------------------------------------------"
    529     END IF
    530     CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
    531   ELSE
    532     CALL grille_m   (lon, lat,        vari, lon1,        lat1, vtmp)
     438  IF(ibeg.AND.prt_level>1) THEN
     439    WRITE(lunout,*)"--------------------------------------------------------"
     440    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
     441    WRITE(lunout,*)"--------------------------------------------------------"
    533442  END IF
     443  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
    534444  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
    535445
Note: See TracChangeset for help on using the changeset viewer.