--- trunk/libf/dyn3d/start_init_orog_m.f90 2008/08/01 15:24:12 15 +++ trunk/Sources/phylmd/Orography/start_init_orog_m.f 2017/10/16 12:35:41 225 @@ -3,116 +3,100 @@ ! From startvar.F, version 1.4 ! 2006/01/27 15:14:22 Fairhead + use dimens_m, only: iim, jjm + IMPLICIT NONE - REAL, ALLOCATABLE, SAVE:: mask(:, :) ! fraction of land (iim + 1, jjm + 1) - REAL, ALLOCATABLE, SAVE:: phis(:, :) ! surface geopotential, in m2 s-2 + REAL, SAVE:: mask(iim + 1, jjm + 1) ! interpolated fraction of land + + private iim, jjm CONTAINS - SUBROUTINE start_init_orog(relief, zstd_2d, zsig_2d, zgam_2d, zthe_2d, & - zpic_2d, zval_2d) + SUBROUTINE start_init_orog(phis, zmea_2d, zstd_2d, zsig_2d, zgam_2d, & + zthe_2d, zpic_2d, zval_2d) - USE ioipsl, only: flininfo, flinopen_nozoom, flinget, flinclo use conf_dat2d_m, only: conf_dat2d - use comgeom, only: rlatu, rlonv - use dimens_m, only: iim, jjm - use indicesol, only: epsfra - use comconst, only: pi + use dynetat0_m, only: rlatu, rlonv use grid_noro_m, only: grid_noro + use indicesol, only: epsfra + use netcdf, only: nf90_nowrite + use netcdf95, only: nf95_open, nf95_gw_var, nf95_inq_varid, nf95_close + use nr_util, only: pi, assert - REAL, intent(out):: relief(:, :) ! orographie moyenne + REAL, intent(out):: phis(:, :) ! (iim + 1, jjm + 1) + ! surface geopotential, in m2 s-2 - REAL, intent(out):: zstd_2d(:, :) + REAL, intent(out):: zmea_2d(:, :) ! (iim + 1, jjm + 1) orographie moyenne + + REAL, intent(out):: zstd_2d(:, :) ! (iim + 1, jjm + 1) ! (deviation standard de l'orographie sous-maille) - REAL, intent(out):: zsig_2d(:, :) + REAL, intent(out):: zsig_2d(:, :) ! (iim + 1, jjm + 1) ! (pente de l'orographie sous-maille) - REAL, intent(out):: zgam_2d(:, :) + REAL, intent(out):: zgam_2d(:, :) ! (iim + 1, jjm + 1) ! (anisotropie de l'orographie sous maille) - REAL, intent(out):: zthe_2d(:, :) + REAL, intent(out):: zthe_2d(:, :) ! (iim + 1, jjm + 1) ! (orientation de l'axe oriente dans la direction de plus grande ! pente de l'orographie sous maille) - REAL, intent(out):: zpic_2d(:, :) ! hauteur pics de la SSO - REAL, intent(out):: zval_2d(:, :) ! hauteur vallees de la SSO + REAL, intent(out):: zpic_2d(:, :) ! (iim + 1, jjm + 1) + ! hauteur pics de la SSO + + REAL, intent(out):: zval_2d(:, :) ! (iim + 1, jjm + 1) + ! hauteur vallees de la SSO ! Local: - INTEGER, SAVE:: iml_rel - INTEGER, SAVE:: jml_rel - REAL lev(1), date, dt - INTEGER itau(1), fid - INTEGER llm_tmp, ttm_tmp - REAL, ALLOCATABLE:: relief_hi(:, :) + INTEGER iml_rel + INTEGER jml_rel + INTEGER ncid, varid + REAL, ALLOCATABLE:: relief(:, :) REAL, ALLOCATABLE:: lon_rad(:), lat_rad(:) REAL, ALLOCATABLE:: lon_ini(:), lat_ini(:) - REAL, ALLOCATABLE:: lon_rel(:, :), lat_rel(:, :) - - CHARACTER(len=120) orogfname !----------------------------------- print *, "Call sequence information: start_init_orog" - if (any((/size(relief, 1), size(zstd_2d, 1), size(zsig_2d, 1), & - size(zgam_2d, 1), size(zthe_2d, 1), size(zpic_2d, 1), & - size(zval_2d, 1)/) /= iim + 1)) stop "start_init_orog size 1" - if (any((/size(relief, 2), size(zstd_2d, 2), size(zsig_2d, 2), & - size(zgam_2d, 2), size(zthe_2d, 2), size(zpic_2d, 2), & - size(zval_2d, 2)/) /= jjm + 1)) stop "start_init_orog size 2" - - orogfname = 'Relief.nc' - print *, 'Reading the high resolution orography' - - CALL flininfo(orogfname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) - - ALLOCATE(lat_rel(iml_rel, jml_rel)) - ALLOCATE(lon_rel(iml_rel, jml_rel)) - ALLOCATE(relief_hi(iml_rel, jml_rel)) - - CALL flinopen_nozoom(orogfname, iml_rel, jml_rel, llm_tmp, & - lon_rel, lat_rel, lev, ttm_tmp, itau, date, dt, fid) - ! 'RELIEF': high resolution orography - CALL flinget(fid, 'RELIEF', iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, & - relief_hi) - CALL flinclo(fid) + call assert((/size(phis, 1), size(zmea_2d, 1), size(zstd_2d, 1), & + size(zsig_2d, 1), size(zgam_2d, 1), size(zthe_2d, 1), & + size(zpic_2d, 1), size(zval_2d, 1)/) == iim + 1, "start_init_orog iim") + call assert((/size(phis, 2), size(zmea_2d, 2), size(zstd_2d, 2), & + size(zsig_2d, 2), size(zgam_2d, 2), size(zthe_2d, 2), & + size(zpic_2d, 2), size(zval_2d, 2)/) == jjm + 1, "start_init_orog jjm") + + print *, 'Reading the high resolution orography...' + + call nf95_open('Relief.nc', nf90_nowrite, ncid) + + call nf95_inq_varid(ncid, "longitude", varid) + call nf95_gw_var(ncid, varid, lon_ini) + lon_ini = lon_ini * pi / 180. ! convert to rad + iml_rel = size(lon_ini) + + call nf95_inq_varid(ncid, "latitude", varid) + call nf95_gw_var(ncid, varid, lat_ini) + lat_ini = lat_ini * pi / 180. ! convert to rad + jml_rel = size(lat_ini) - ! In case we have a file which is in degrees we do the transformation: + call nf95_inq_varid(ncid, "RELIEF", varid) + call nf95_gw_var(ncid, varid, relief) - ALLOCATE(lon_rad(iml_rel)) - ALLOCATE(lon_ini(iml_rel)) - - IF (MAXVAL(lon_rel(:, :)) > pi) THEN - lon_ini(:) = lon_rel(:, 1) * pi / 180. - ELSE - lon_ini(:) = lon_rel(:, 1) - ENDIF + call nf95_close(ncid) + ALLOCATE(lon_rad(iml_rel)) ALLOCATE(lat_rad(jml_rel)) - ALLOCATE(lat_ini(jml_rel)) - IF (MAXVAL(lat_rel(:, :)) > pi) THEN - lat_ini(:) = lat_rel(1, :) * pi / 180. - ELSE - lat_ini(:) = lat_rel(1, :) - ENDIF - - CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, relief_hi , & + CALL conf_dat2d(lon_ini, lat_ini, lon_rad, lat_rad, relief , & interbar=.FALSE.) print *, 'Compute all the parameters needed for the gravity wave drag code' - ! Allocate the data we need to put in the interpolated fields: - ALLOCATE(phis(iim + 1, jjm + 1)) - ALLOCATE(mask(iim + 1, jjm + 1)) - - CALL grid_noro(lon_rad, lat_rad, relief_hi, rlonv, rlatu, phis, relief, & + CALL grid_noro(lon_rad, lat_rad, relief, rlonv, rlatu, phis, zmea_2d, & zstd_2d, zsig_2d, zgam_2d, zthe_2d, zpic_2d, zval_2d, mask) - - phis(iim + 1, :) = phis(1, :) phis(:, :) = phis(:, :) * 9.81 mask(2:, 1) = mask(1, 1) ! north pole