module fyhyp_m IMPLICIT NONE contains SUBROUTINE fyhyp(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1) ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32 ! Author: P. Le Van, from analysis by R. Sadourny ! 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 include "dimensions.h" ! for jjm include "serre.h" ! for clat, grossismy, dzoomy, tauy REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1) REAL, intent(out):: rlatv(jjm) real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm) ! Local: DOUBLE PRECISION champmin, champmax 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, y0, pisjm DOUBLE PRECISION yo1, yi, ylon2, ymoy, yprimin DOUBLE PRECISION yfi, yf1, ffdy DOUBLE PRECISION ypn, deply, y00 SAVE y00, deply INTEGER i, j, it, ik, iter, jlat INTEGER jpn, jjpn SAVE jpn DOUBLE PRECISION a0, a1, a2, a3, yi2, heavyy0, heavyy0m DOUBLE PRECISION fa(0:nmax2), fb(0:nmax2) REAL y0min, y0max DOUBLE PRECISION heavyside !------------------------------------------------------------------- print *, "Call sequence information: fyhyp" pi = 2.*asin(1.) pis2 = pi/2. pisjm = pi/real(jjm) epsilon = 1e-3 y0 = clat*pi/180. dzoom = dzoomy*pi print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):' print *, y0, grossismy, tauy, dzoom DO i = 0, nmax2 yt(i) = -pis2 + real(i)*pi/nmax2 END DO heavyy0m = heavyside(-y0) heavyy0 = heavyside(y0) y0min = 2.*y0*heavyy0m - pis2 y0max = 2.*y0*heavyy0 + pis2 fa = 999.999 fb = 999.999 DO i = 0, nmax2 IF (yt(i)y0) THEN fa(i) = tauy*(y0-yt(i) + dzoom/2.) fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0) END IF IF (200.*fb(i)<-fa(i)) THEN fhyp(i) = -1. ELSE IF (200.*fb(i)y0) THEN fa(i) = tauy*(y0-ymoy + dzoom/2.) fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0) 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), a0, a1, a2, a3) 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, champmin, champmax 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 ') END SUBROUTINE fyhyp end module fyhyp_m