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).
File:
1 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!-------------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.