--- trunk/libf/dyn3d/dynetat0.f90 2008/07/25 19:59:34 13 +++ trunk/dyn3d/dynetat0.f 2018/07/12 17:53:18 278 @@ -1,193 +1,858 @@ module dynetat0_m - ! This module is clean: no C preprocessor directive, no include line. + use dimensions, only: iim, jjm IMPLICIT NONE -contains + private iim, jjm, principal_cshift, invert_zoom_x, funcd + + INTEGER day_ini + ! day number at the beginning of the run, based at value 1 on + ! January 1st of annee_ref + + integer:: day_ref = 1 ! jour de l'ann\'ee de l'\'etat initial + ! (= 350 si 20 d\'ecembre par exemple) + + integer:: annee_ref = 1998 ! Annee de l'etat initial (avec 4 chiffres) + + REAL, protected:: clon ! longitude of the center of the zoom, in rad + real, protected:: clat ! latitude of the center of the zoom, in rad + + real, protected:: grossismx, grossismy + ! facteurs de grossissement du zoom, selon la longitude et la latitude + ! = 2 si 2 fois, = 3 si 3 fois, etc. + + real, protected:: dzoomx, dzoomy + ! extensions en longitude et latitude de la zone du zoom (fractions + ! de la zone totale) + + real, protected:: taux, tauy + ! raideur de la transition de l'int\'erieur \`a l'ext\'erieur du zoom + + real, protected:: rlatu(jjm + 1) + ! latitudes of points of the "scalar" and "u" grid, in rad + + real, protected:: rlatv(jjm) + ! latitudes of points of the "v" grid, in rad, in decreasing order + + real, protected:: rlonu(iim + 1) + ! longitudes of points of the "u" grid, in rad - SUBROUTINE dynetat0(vcov, ucov, teta, q, masse, ps, phis, time) + real, protected:: rlonv(iim + 1) + ! longitudes of points of the "scalar" and "v" grid, in rad - ! From dynetat0.F, version 1.2 2004/06/22 11:45:30 + real, protected:: xprimu(iim + 1), xprimv(iim + 1) + ! 2 pi / iim * (derivative of the longitudinal zoom function)(rlon[uv]) - ! Authors: P. Le Van, L. Fairhead - ! Objet : lecture de l'état initial + REAL, protected:: xprimm025(iim + 1), xprimp025(iim + 1) + REAL, protected:: rlatu1(jjm), rlatu2(jjm), yprimu1(jjm), yprimu2(jjm) + REAL ang0, etot0, ptot0, ztot0, stot0 + INTEGER, PARAMETER, private:: nmax = 30000 + DOUBLE PRECISION, private:: abs_y - use dimens_m, only: iim, jjm, llm, nqmx - use comconst, only: im, cpp, dtvr, g, kappa, jm, lllm, omeg, rad - use comvert, only: pa - use logic, only: fxyhypb, ysinus - use comgeom, only: rlonu, rlatu, rlonv, rlatv, cu_2d, cv_2d, aire_2d - use serre, only: clon, clat, grossismy, grossismx - use temps, only: day_ref, day_ini, itau_dyn, annee_ref - use ener, only: etot0, ang0, ptot0, stot0, ztot0 - use advtrac_m, only: tname - use netcdf95, only: nf95_open, nf95_inq_varid, handle_err, NF95_CLOSE - use netcdf, only: NF90_NOWRITE, NF90_GET_VAR, NF90_NOERR - use numer_rec, only: assert - - ! Arguments: - REAL, intent(out):: vcov(: , :), ucov(:, :), teta(:, :) - REAL, intent(out):: q(:, :, :), masse(:, :) - REAL, intent(out):: ps(:) ! in Pa - REAL, intent(out):: phis(:, :) - REAL, intent(out):: time - - ! Variables - INTEGER length, iq - PARAMETER (length = 100) - REAL tab_cntrl(length) ! tableau des parametres du run - INTEGER ierr, nid, nvarid + save + +contains + + SUBROUTINE dynetat0(vcov, ucov, teta, q, masse, ps, phis) + + ! From dynetat0.F, version 1.2, 2004/06/22 11:45:30 + ! Authors: P. Le Van, L. Fairhead + ! This procedure reads the initial state of the atmosphere. + + use comconst, only: dtvr + use conf_gcm_m, only: raz_date + use dimensions, only: iim, jjm, llm, nqmx + use disvert_m, only: pa + use iniadvtrac_m, only: tname + use netcdf, only: NF90_NOWRITE, NF90_NOERR + use netcdf95, only: NF95_GET_VAR, nf95_open, nf95_inq_varid, NF95_CLOSE, & + NF95_Gw_VAR + use nr_util, only: assert + use temps, only: itau_dyn + use unit_nml_m, only: unit_nml + + REAL, intent(out):: vcov(: , :, :) ! (iim + 1, jjm, llm) + REAL, intent(out):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) + REAL, intent(out):: teta(:, :, :) ! (iim + 1, jjm + 1, llm) + REAL, intent(out):: q(:, :, :, :) ! (iim + 1, jjm + 1, llm, nqmx) + REAL, intent(out):: masse(:, :, :) ! (iim + 1, jjm + 1, llm) + REAL, intent(out):: ps(:, :) ! (iim + 1, jjm + 1) in Pa + REAL, intent(out):: phis(:, :) ! (iim + 1, jjm + 1) + + ! Local variables: + INTEGER iq + REAL, allocatable:: tab_cntrl(:) ! tableau des param\`etres du run + INTEGER ierr, ncid, varid + + namelist /dynetat0_nml/ day_ref, annee_ref !----------------------------------------------------------------------- print *, "Call sequence information: dynetat0" - call assert(size(vcov, 1) == (iim + 1) * jjm, "dynetat0 vcov 1") - call assert((/size(ucov, 1), size(teta, 1), size(q, 1), size(masse, 1), & - size(ps)/) == (iim + 1) * (jjm + 1), "dynetat0 (iim + 1) * (jjm + 1)") - call assert(shape(phis) == (/iim + 1, jjm + 1/), "dynetat0 phis") - call assert((/size(vcov, 2), size(ucov, 2), size(teta, 2), size(q, 2), & - size(masse, 2)/) == llm, "dynetat0 llm") - call assert(size(q, 3) == nqmx, "dynetat0 q 3") - - ! Fichier état initial : - call nf95_open("start.nc", NF90_NOWRITE, nid) - - call nf95_inq_varid(nid, "controle", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, tab_cntrl) - call handle_err("dynetat0, controle", ierr, nid) - - im = int(tab_cntrl(1)) - jm = int(tab_cntrl(2)) - lllm = int(tab_cntrl(3)) - day_ref = int(tab_cntrl(4)) - annee_ref = int(tab_cntrl(5)) - omeg = tab_cntrl(7) - dtvr = tab_cntrl(12) - etot0 = tab_cntrl(13) - ptot0 = tab_cntrl(14) - ztot0 = tab_cntrl(15) - stot0 = tab_cntrl(16) - ang0 = tab_cntrl(17) - pa = tab_cntrl(18) - clon = tab_cntrl(20) - clat = tab_cntrl(21) - grossismx = tab_cntrl(22) - grossismy = tab_cntrl(23) - - IF (tab_cntrl(24) == 1.) THEN - fxyhypb = .TRUE. - ELSE - fxyhypb = .FALSE. - ysinus = .FALSE. - IF (tab_cntrl(27) == 1.) ysinus = .TRUE. + call assert((/size(ucov, 1), size(vcov, 1), size(masse, 1), size(ps, 1), & + size(phis, 1), size(q, 1), size(teta, 1)/) == iim + 1, "dynetat0 iim") + call assert((/size(ucov, 2), size(vcov, 2) + 1, size(masse, 2), & + size(ps, 2), size(phis, 2), size(q, 2), size(teta, 2)/) == jjm + 1, & + "dynetat0 jjm") + call assert((/size(vcov, 3), size(ucov, 3), size(teta, 3), size(q, 3), & + size(masse, 3)/) == llm, "dynetat0 llm") + call assert(size(q, 4) == nqmx, "dynetat0 q nqmx") + + ! Fichier \'etat initial : + call nf95_open("start.nc", NF90_NOWRITE, ncid) + + call nf95_inq_varid(ncid, "controle", varid) + call NF95_Gw_VAR(ncid, varid, tab_cntrl) + + call assert(int(tab_cntrl(1)) == iim, "dynetat0 tab_cntrl iim") + call assert(int(tab_cntrl(2)) == jjm, "dynetat0 tab_cntrl jjm") + call assert(int(tab_cntrl(3)) == llm, "dynetat0 tab_cntrl llm") + + IF (dtvr /= tab_cntrl(12)) THEN + print *, 'Warning: the time steps from day_step and "start.nc" ' // & + 'are different.' + print *, 'dtvr from day_step: ', dtvr + print *, 'dtvr from "start.nc": ', tab_cntrl(12) + print *, 'Using the value from day_step.' ENDIF - day_ini = tab_cntrl(30) - itau_dyn = tab_cntrl(31) + etot0 = tab_cntrl(13) + ptot0 = tab_cntrl(14) + ztot0 = tab_cntrl(15) + stot0 = tab_cntrl(16) + ang0 = tab_cntrl(17) + pa = tab_cntrl(18) + + clon = tab_cntrl(20) + clat = tab_cntrl(21) + grossismx = tab_cntrl(22) + grossismy = tab_cntrl(23) + dzoomx = tab_cntrl(25) + dzoomy = tab_cntrl(26) + taux = tab_cntrl(28) + tauy = tab_cntrl(29) + + print *, "Enter namelist 'dynetat0_nml'." + read(unit=*, nml=dynetat0_nml) + write(unit_nml, nml=dynetat0_nml) + + if (raz_date) then + print *, 'Resetting the date, using the namelist.' + day_ini = day_ref + itau_dyn = 0 + else + day_ref = tab_cntrl(4) + annee_ref = tab_cntrl(5) + itau_dyn = tab_cntrl(31) + day_ini = tab_cntrl(30) + end if + + print *, "day_ini = ", day_ini + + call NF95_INQ_VARID (ncid, "rlonu", varid) + call NF95_GET_VAR(ncid, varid, rlonu) + + call NF95_INQ_VARID (ncid, "rlatu", varid) + call NF95_GET_VAR(ncid, varid, rlatu) + + call NF95_INQ_VARID (ncid, "rlonv", varid) + call NF95_GET_VAR(ncid, varid, rlonv) + + call NF95_INQ_VARID (ncid, "rlatv", varid) + call NF95_GET_VAR(ncid, varid, rlatv) + + CALL nf95_inq_varid(ncid, 'xprimu', varid) + CALL nf95_get_var(ncid, varid, xprimu) + + CALL nf95_inq_varid(ncid, 'xprimv', varid) + CALL nf95_get_var(ncid, varid, xprimv) + + CALL nf95_inq_varid(ncid, 'xprimm025', varid) + CALL nf95_get_var(ncid, varid, xprimm025) + + CALL nf95_inq_varid(ncid, 'xprimp025', varid) + CALL nf95_get_var(ncid, varid, xprimp025) + + call NF95_INQ_VARID (ncid, "rlatu1", varid) + call NF95_GET_VAR(ncid, varid, rlatu1) + + call NF95_INQ_VARID (ncid, "rlatu2", varid) + call NF95_GET_VAR(ncid, varid, rlatu2) + + CALL nf95_inq_varid(ncid, 'yprimu1', varid) + CALL nf95_get_var(ncid, varid, yprimu1) + + CALL nf95_inq_varid(ncid, 'yprimu2', varid) + CALL nf95_get_var(ncid, varid, yprimu2) + + call NF95_INQ_VARID (ncid, "phis", varid) + call NF95_GET_VAR(ncid, varid, phis) - PRINT *, 'rad = ', rad - PRINT *, 'omeg = ', omeg - PRINT *, 'g = ', g - PRINT *, 'cpp = ', cpp - PRINT *, 'kappa = ', kappa + call NF95_INQ_VARID (ncid, "ucov", varid) + call NF95_GET_VAR(ncid, varid, ucov) - IF (im /= iim) THEN - PRINT 1, im, iim - STOP 1 - ELSE IF (jm /= jjm) THEN - PRINT 2, jm, jjm - STOP 1 - ELSE IF (lllm /= llm) THEN - PRINT 3, lllm, llm - STOP 1 - ENDIF + call NF95_INQ_VARID (ncid, "vcov", varid) + call NF95_GET_VAR(ncid, varid, vcov) - call NF95_INQ_VARID (nid, "rlonu", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, rlonu) - call handle_err("dynetat0, rlonu", ierr, nid) - - call NF95_INQ_VARID (nid, "rlatu", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, rlatu) - call handle_err("dynetat0, rlatu", ierr, nid) - - call NF95_INQ_VARID (nid, "rlonv", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, rlonv) - call handle_err("dynetat0, rlonv", ierr, nid) - - call NF95_INQ_VARID (nid, "rlatv", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, rlatv) - call handle_err("dynetat0, rlatv", ierr, nid) - - call NF95_INQ_VARID (nid, "cu", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, cu_2d) - call handle_err("dynetat0, cu", ierr, nid) - - call NF95_INQ_VARID (nid, "cv", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, cv_2d) - call handle_err("dynetat0, cv", ierr, nid) - - call NF95_INQ_VARID (nid, "aire", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, aire_2d) - call handle_err("dynetat0, aire", ierr, nid) - - call NF95_INQ_VARID (nid, "phisinit", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, phis) - call handle_err("dynetat0, phisinit", ierr, nid) - - call NF95_INQ_VARID (nid, "temps", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, time) - call handle_err("dynetat0, temps", ierr, nid) - - call NF95_INQ_VARID (nid, "ucov", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, ucov, count=(/iim + 1, jjm + 1, llm/)) - call handle_err("dynetat0, ucov", ierr, nid) - - call NF95_INQ_VARID (nid, "vcov", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, vcov, count=(/iim + 1, jjm, llm/)) - call handle_err("dynetat0, vcov", ierr, nid) - - call NF95_INQ_VARID (nid, "teta", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, teta, count=(/iim + 1, jjm + 1, llm/)) - call handle_err("dynetat0, teta", ierr, nid) + call NF95_INQ_VARID (ncid, "teta", varid) + call NF95_GET_VAR(ncid, varid, teta) DO iq = 1, nqmx - call NF95_INQ_VARID(nid, tname(iq), nvarid, ierr) - IF (ierr /= NF90_NOERR) THEN - PRINT *, 'dynetat0: le champ "' // tname(iq) // '" est absent, ' // & - "il est donc initialisé à zéro." - q(:, :, iq) = 0. + call NF95_INQ_VARID(ncid, tname(iq), varid, ierr) + IF (ierr == NF90_NOERR) THEN + call NF95_GET_VAR(ncid, varid, q(:, :, :, iq)) ELSE - ierr = NF90_GET_VAR(nid, nvarid, q(:, :, iq), & - count=(/iim + 1, jjm + 1, llm/)) - call handle_err("dynetat0, " // tname(iq), ierr, nid) + PRINT *, 'dynetat0: "' // tname(iq) // '" not found, ' // & + "setting it to zero..." + q(:, :, :, iq) = 0. ENDIF ENDDO - call NF95_INQ_VARID (nid, "masse", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, masse, count=(/iim + 1, jjm + 1, llm/)) - call handle_err("dynetat0, masse", ierr, nid) - - call NF95_INQ_VARID (nid, "ps", nvarid) - ierr = NF90_GET_VAR(nid, nvarid, ps, count=(/iim + 1, jjm + 1/)) - call handle_err("dynetat0, ps", ierr, nid) - - call NF95_CLOSE(nid) - - day_ini=day_ini+INT(time) - time=time-INT(time) - -1 FORMAT(//10x, 'la valeur de im =', i4, 2x, & - 'lue sur le fichier de demarrage est differente de la valeur ' & - // 'parametree iim =', i4//) -2 FORMAT(//10x, 'la valeur de jm =', i4, 2x, & - 'lue sur le fichier de demarrage est differente de la valeur ' & - // 'parametree jjm =', i4//) -3 FORMAT(//10x, 'la valeur de lmax =', i4, 2x, & - 'lue sur le fichier demarrage est differente de la valeur ' & - // 'parametree llm =', i4//) + call NF95_INQ_VARID (ncid, "masse", varid) + call NF95_GET_VAR(ncid, varid, masse) + + call NF95_INQ_VARID (ncid, "ps", varid) + call NF95_GET_VAR(ncid, varid, ps) + ! Check that there is a single value at each pole: + call assert(ps(1, 1) == ps(2:, 1), "dynetat0 ps north pole") + call assert(ps(1, jjm + 1) == ps(2:, jjm + 1), "dynetat0 ps south pole") + + call NF95_CLOSE(ncid) END SUBROUTINE dynetat0 + !******************************************************************** + + subroutine read_serre + + use unit_nml_m, only: unit_nml + use nr_util, only: assert, pi + + REAL:: clon_deg = 0. ! longitude of the center of the zoom, in degrees + real:: clat_deg = 0. ! latitude of the center of the zoom, in degrees + + namelist /serre_nml/ clon_deg, clat_deg, grossismx, grossismy, dzoomx, & + dzoomy, taux, tauy + + !------------------------------------------------- + + ! Default values: + grossismx = 1. + grossismy = 1. + dzoomx = 0.2 + dzoomy = 0.2 + taux = 3. + tauy = 3. + + print *, "Enter namelist 'serre_nml'." + read(unit=*, nml=serre_nml) + write(unit_nml, nml=serre_nml) + + call assert(grossismx >= 1. .and. grossismy >= 1., "read_serre grossism") + call assert(dzoomx > 0., dzoomx < 1., dzoomy < 1., & + "read_serre dzoomx dzoomy") + clon = clon_deg / 180. * pi + clat = clat_deg / 180. * pi + + end subroutine read_serre + + !******************************************************************** + + SUBROUTINE fyhyp + + ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32 + + ! Author: P. Le Van, from analysis by R. Sadourny + + ! Define rlatu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1, using + ! clat, grossismy, dzoomy, tauy. + + ! Calcule les latitudes et dérivées dans la grille du GCM pour une + ! fonction f(y) à dérivée tangente hyperbolique. + + ! Il vaut mieux avoir : grossismy * dzoom < pi / 2 + + use coefpoly_m, only: coefpoly, a0, a1, a2, a3 + USE dimensions, only: jjm + use heavyside_m, only: heavyside + + ! Local: + + INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax + REAL dzoom ! distance totale de la zone du zoom (en radians) + DOUBLE PRECISION ylat(jjm + 1), yprim(jjm + 1) + DOUBLE PRECISION yuv + DOUBLE PRECISION, save:: yt(0:nmax2) + DOUBLE PRECISION fhyp(0:nmax2), beta + DOUBLE PRECISION, save:: ytprim(0:nmax2) + DOUBLE PRECISION fxm(0:nmax2) + DOUBLE PRECISION, save:: yf(0:nmax2) + DOUBLE PRECISION yypr(0:nmax2) + DOUBLE PRECISION yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1) + DOUBLE PRECISION pi, pis2, epsilon, pisjm + DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin + DOUBLE PRECISION yfi, yf1, ffdy + DOUBLE PRECISION ypn + DOUBLE PRECISION, save::deply, y00 + + INTEGER i, j, it, ik, iter, jlat, jjpn + INTEGER, save:: jpn + DOUBLE PRECISION yi2, heavyy0, heavyy0m + DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2) + REAL y0min, y0max + + !------------------------------------------------------------------- + + print *, "Call sequence information: fyhyp" + + pi = 2.*asin(1.) + pis2 = pi/2. + pisjm = pi/real(jjm) + epsilon = 1e-3 + dzoom = dzoomy*pi + + DO i = 0, nmax2 + yt(i) = -pis2 + real(i)*pi/nmax2 + END DO + + heavyy0m = heavyside(-clat) + heavyy0 = heavyside(clat) + y0min = 2.*clat*heavyy0m - pis2 + y0max = 2.*clat*heavyy0 + pis2 + + fa = 999.999 + fb = 999.999 + + DO i = 0, nmax2 + IF (yt(i)clat) THEN + fa(i) = tauy*(clat-yt(i) + dzoom/2.) + fb(i) = (2.*clat*heavyy0-yt(i) + pis2)*(yt(i)-clat) + END IF + + IF (200.*fb(i)<-fa(i)) THEN + fhyp(i) = -1. + ELSE IF (200.*fb(i)clat) THEN + fa(i) = tauy*(clat-ymoy + dzoom/2.) + fb(i) = (2.*clat*heavyy0-ymoy + pis2)*(ymoy-clat) + END IF + + IF (200.*fb(i)<-fa(i)) THEN + fxm(i) = -1. + ELSE IF (200.*fb(i)= 1 .and. yfi < yf(it)) + it = it - 1 + END DO + + yi = yt(it) + IF (it==nmax2) THEN + it = nmax2 - 1 + yf(it + 1) = pis2 + END IF + + ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi) + ! et Y'(yi) + + CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), & + yt(it), yt(it + 1)) + + yf1 = yf(it) + yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi + + iter = 1 + DO + yi = yi - (yf1-yfi)/yprimin + IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit + yo1 = yi + yi2 = yi*yi + yf1 = a0 + a1*yi + a2*yi2 + a3*yi2*yi + yprimin = a1 + 2.*a2*yi + 3.*a3*yi2 + END DO + if (abs(yi-yo1) > epsilon) then + print *, 'Pas de solution.', j, ylon2 + STOP 1 + end if + + yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi + yprim(j) = pi/(jjm*yprimin) + yvrai(j) = yi + END DO + + DO j = 1, jlat - 1 + IF (yvrai(j + 1)= rlatu1(j)) THEN + print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j + STOP 17 + ENDIF + + IF (rlatv(j) >= rlatu(j)) THEN + print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j + STOP 18 + ENDIF + ENDDO + + print *, 'Latitudes' + print 3, minval(ylat(:jjm)) *180d0/pi, maxval(ylat(:jjm))*180d0/pi + +3 Format(1x, ' Au centre du zoom, la longueur de la maille est', & + ' d environ ', f0.2, ' degres ', /, & + ' alors que la maille en dehors de la zone du zoom est ', & + "d'environ ", f0.2, ' degres ') + + rlatu(1) = pi / 2. + rlatu(jjm + 1) = -rlatu(1) + + END SUBROUTINE fyhyp + + !******************************************************************** + + SUBROUTINE fxhyp + + ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32 + ! Author: P. Le Van, from formulas by R. Sadourny + + ! Compute xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025, + ! using clon, grossismx, dzoomx, taux. + + ! Calcule les longitudes et dérivées dans la grille du GCM pour + ! une fonction x_f(\tilde x) à dérivée tangente hyperbolique. + + ! Il vaut mieux avoir : grossismx \times delta < pi + + ! Le premier point scalaire pour une grille regulière (grossismx = + ! 1) avec clon = 0 est à - 180 degrés. + + use nr_util, only: pi, pi_d, twopi, twopi_d, arth, assert, rad_to_deg + + USE dimensions, ONLY: iim + use tanh_cautious_m, only: tanh_cautious + + ! Local: + real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim) + REAL delta, h + DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf, ffdx + DOUBLE PRECISION beta + INTEGER i, is2 + DOUBLE PRECISION xmoy(nmax), fxm(nmax) + + !---------------------------------------------------------------------- + + print *, "Call sequence information: fxhyp" + + if (grossismx == 1.) then + h = twopi / iim + + xprimm025(:iim) = h + xprimp025(:iim) = h + xprimv(:iim) = h + xprimu(:iim) = h + + rlonv(:iim) = arth(- pi + clon, h, iim) + rlonm025(:iim) = rlonv(:iim) - 0.25 * h + rlonp025(:iim) = rlonv(:iim) + 0.25 * h + rlonu(:iim) = rlonv(:iim) + 0.5 * h + else + delta = dzoomx * twopi_d + xtild = arth(0d0, pi_d / nmax, nmax + 1) + forall (i = 1:nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i)) + + ! Compute fhyp: + fhyp(1:nmax - 1) = tanh_cautious(taux * (delta / 2d0 & + - xtild(1:nmax - 1)), xtild(1:nmax - 1) & + * (pi_d - xtild(1:nmax - 1))) + fhyp(0) = 1d0 + fhyp(nmax) = -1d0 + + fxm = tanh_cautious(taux * (delta / 2d0 - xmoy), xmoy * (pi_d - xmoy)) + + ! Compute \int_0 ^{\tilde x} F: + + ffdx(0) = 0d0 + + DO i = 1, nmax + ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1)) + END DO + + print *, "ffdx(nmax) = ", ffdx(nmax) + beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax)) + print *, "beta = ", beta + + IF (2d0 * beta - grossismx <= 0d0) THEN + print *, 'Bad choice of grossismx, taux, dzoomx.' + print *, 'Decrease dzoomx or grossismx.' + STOP 1 + END IF + + G = beta + (grossismx - beta) * fhyp + + Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) & + * ffdx(:nmax - 1) + Xf(nmax) = pi_d + + call invert_zoom_x(beta, xf, xtild, G, rlonm025(:iim), xprimm025(:iim), & + xuv = - 0.25d0) + call invert_zoom_x(beta, xf, xtild, G, rlonv(:iim), xprimv(:iim), & + xuv = 0d0) + call invert_zoom_x(beta, xf, xtild, G, rlonu(:iim), xprimu(:iim), & + xuv = 0.5d0) + call invert_zoom_x(beta, xf, xtild, G, rlonp025(:iim), xprimp025(:iim), & + xuv = 0.25d0) + end if + + is2 = 0 + + IF (MINval(rlonm025(:iim)) < - pi - 0.1 & + .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN + IF (clon <= 0.) THEN + is2 = 1 + + do while (rlonm025(is2) < - pi .and. is2 < iim) + is2 = is2 + 1 + end do + + call assert(rlonm025(is2) >= - pi, & + "fxhyp -- rlonm025 should be >= - pi") + ELSE + is2 = iim + + do while (rlonm025(is2) > pi .and. is2 > 1) + is2 = is2 - 1 + end do + + if (rlonm025(is2) > pi) then + print *, 'Rlonm025 plus grand que pi !' + STOP 1 + end if + END IF + END IF + + call principal_cshift(is2, rlonm025, xprimm025) + call principal_cshift(is2, rlonv, xprimv) + call principal_cshift(is2, rlonu, xprimu) + call principal_cshift(is2, rlonp025, xprimp025) + + forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i) + print *, "Minimum longitude step:", MINval(d_rlonv) * rad_to_deg, "degrees" + print *, "Maximum longitude step:", MAXval(d_rlonv) * rad_to_deg, "degrees" + + ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu: + DO i = 1, iim + 1 + IF (rlonp025(i) < rlonv(i)) THEN + print *, 'rlonp025(', i, ') = ', rlonp025(i) + print *, "< rlonv(", i, ") = ", rlonv(i) + STOP 1 + END IF + + IF (rlonv(i) < rlonm025(i)) THEN + print *, 'rlonv(', i, ') = ', rlonv(i) + print *, "< rlonm025(", i, ") = ", rlonm025(i) + STOP 1 + END IF + + IF (rlonp025(i) > rlonu(i)) THEN + print *, 'rlonp025(', i, ') = ', rlonp025(i) + print *, "> rlonu(", i, ") = ", rlonu(i) + STOP 1 + END IF + END DO + + END SUBROUTINE fxhyp + + !******************************************************************** + + subroutine principal_cshift(is2, xlon, xprimm) + + ! Add or subtract 2 pi so that xlon is near [-pi, pi], then cshift + ! so that xlon is in ascending order. Make the same cshift on + ! xprimm. Use clon. + + use nr_util, only: twopi + + USE dimensions, ONLY: iim + + integer, intent(in):: is2 + real, intent(inout):: xlon(:), xprimm(:) ! (iim + 1) + + !----------------------------------------------------- + + if (is2 /= 0) then + IF (clon <= 0.) THEN + IF (is2 /= 1) THEN + xlon(:is2 - 1) = xlon(:is2 - 1) + twopi + xlon(:iim) = cshift(xlon(:iim), shift = is2 - 1) + xprimm(:iim) = cshift(xprimm(:iim), shift = is2 - 1) + END IF + else + xlon(is2 + 1:iim) = xlon(is2 + 1:iim) - twopi + xlon(:iim) = cshift(xlon(:iim), shift = is2) + xprimm(:iim) = cshift(xprimm(:iim), shift = is2) + end IF + end if + + xlon(iim + 1) = xlon(1) + twopi + xprimm(iim + 1) = xprimm(1) + + end subroutine principal_cshift + + !********************************************************************** + + subroutine invert_zoom_x(beta, xf, xtild, G, xlon, xprim, xuv) + + ! Using clon and grossismx. + + use coefpoly_m, only: coefpoly, a1, a2, a3 + USE dimensions, ONLY: iim + use nr_util, only: pi_d, twopi_d + use numer_rec_95, only: hunt, rtsafe + + DOUBLE PRECISION, intent(in):: beta, Xf(0:), xtild(0:), G(0:) ! (0:nmax) + + real, intent(out):: xlon(:), xprim(:) ! (iim) + + DOUBLE PRECISION, intent(in):: xuv + ! between - 0.25 and 0.5 + ! 0. si calcul aux points scalaires + ! 0.5 si calcul aux points U + + ! Local: + DOUBLE PRECISION Y + DOUBLE PRECISION h ! step of the uniform grid + integer i, it + + DOUBLE PRECISION xvrai(iim), Gvrai(iim) + ! intermediary variables because xlon and xprim are single precision + + !------------------------------------------------------------------ + + print *, "Call sequence information: invert_zoom_x" + it = 0 ! initial guess + h = twopi_d / iim + + DO i = 1, iim + Y = - pi_d + (i + xuv - 0.75d0) * h + ! - pi <= y < pi + abs_y = abs(y) + + ! Distinguish boundaries in order to avoid roundoff error. + ! funcd should be exactly equal to 0 at xtild(it) or xtild(it + + ! 1) and could be very small with the wrong sign so rtsafe + ! would fail. + if (abs_y == 0d0) then + xvrai(i) = 0d0 + gvrai(i) = grossismx + else if (abs_y == pi_d) then + xvrai(i) = pi_d + gvrai(i) = 2d0 * beta - grossismx + else + call hunt(xf, abs_y, it, my_lbound = 0) + ! {0 <= it <= nmax - 1} + + ! Calcul de xvrai(i) et Gvrai(i) + CALL coefpoly(Xf(it), Xf(it + 1), G(it), G(it + 1), xtild(it), & + xtild(it + 1)) + xvrai(i) = rtsafe(funcd, xtild(it), xtild(it + 1), xacc = 1d-6) + Gvrai(i) = a1 + xvrai(i) * (2d0 * a2 + xvrai(i) * 3d0 * a3) + end if + + if (y < 0d0) xvrai(i) = - xvrai(i) + end DO + + DO i = 1, iim -1 + IF (xvrai(i + 1) < xvrai(i)) THEN + print *, 'xvrai(', i + 1, ') < xvrai(', i, ')' + STOP 1 + END IF + END DO + + xlon = xvrai + clon + xprim = h / Gvrai + + end subroutine invert_zoom_x + + !********************************************************************** + + SUBROUTINE funcd(x, fval, fderiv) + + use coefpoly_m, only: a0, a1, a2, a3 + + DOUBLE PRECISION, INTENT(IN):: x + DOUBLE PRECISION, INTENT(OUT):: fval, fderiv + + fval = a0 + x * (a1 + x * (a2 + x * a3)) - abs_y + fderiv = a1 + x * (2d0 * a2 + x * 3d0 * a3) + + END SUBROUTINE funcd + end module dynetat0_m