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/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.