Index: LMDZ5/branches/testing/libf/dyn3d_common/coefpoly.F =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/coefpoly.F (revision 2220) +++ (revision ) @@ -1,40 +1,0 @@ -! -! $Header$ -! - SUBROUTINE coefpoly ( Xf1, Xf2, Xprim1, Xprim2, xtild1,xtild2 , - , a0,a1,a2,a3 ) - IMPLICIT NONE -c -c ... Auteur : P. Le Van ... -c -c -c Calcul des coefficients a0, a1, a2, a3 du polynome de degre 3 qui -c satisfait aux 4 equations suivantes : - -c a0 + a1*xtild1 + a2*xtild1*xtild1 + a3*xtild1*xtild1*xtild1 = Xf1 -c a0 + a1*xtild2 + a2*xtild2*xtild2 + a3*xtild2*xtild2*xtild2 = Xf2 -c a1 + 2.*a2*xtild1 + 3.*a3*xtild1*xtild1 = Xprim1 -c a1 + 2.*a2*xtild2 + 3.*a3*xtild2*xtild2 = Xprim2 - -c On en revient a resoudre un systeme de 4 equat.a 4 inconnues a0,a1,a2,a3 - - REAL(KIND=8) Xf1, Xf2,Xprim1,Xprim2, xtild1,xtild2, xi - REAL(KIND=8) Xfout, Xprim - REAL(KIND=8) a1,a2,a3,a0, xtil1car, xtil2car,derr,x1x2car - - xtil1car = xtild1 * xtild1 - xtil2car = xtild2 * xtild2 - - derr= 2. *(Xf2-Xf1)/( xtild1-xtild2) - - x1x2car = ( xtild1-xtild2)*(xtild1-xtild2) - - a3 = (derr + Xprim1+Xprim2 )/x1x2car - a2 = ( Xprim1 - Xprim2 + 3.* a3 * ( xtil2car-xtil1car ) ) / - / ( 2.* ( xtild1 - xtild2 ) ) - - a1 = Xprim1 -3.* a3 * xtil1car -2.* a2 * xtild1 - a0 = Xf1 - a3 * xtild1* xtil1car -a2 * xtil1car - a1 *xtild1 - - RETURN - END Index: LMDZ5/branches/testing/libf/dyn3d_common/coefpoly_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/coefpoly_m.F90 (revision 2258) +++ LMDZ5/branches/testing/libf/dyn3d_common/coefpoly_m.F90 (revision 2258) @@ -0,0 +1,52 @@ +module coefpoly_m + + IMPLICIT NONE + +contains + + SUBROUTINE coefpoly(xf1, xf2, xprim1, xprim2, xtild1, xtild2, a0, a1, a2, a3) + + ! From LMDZ4/libf/dyn3d/coefpoly.F, version 1.1.1.1 2004/05/19 12:53:05 + + ! Author: P. Le Van + + ! Calcul des coefficients a0, a1, a2, a3 du polynôme de degré 3 qui + ! satisfait aux 4 équations suivantes : + + ! a0 + a1 * xtild1 + a2 * xtild1**2 + a3 * xtild1**3 = Xf1 + ! a0 + a1 * xtild2 + a2 * xtild2**2 + a3 * xtild2**3 = Xf2 + ! a1 + 2. * a2 * xtild1 + 3. * a3 * xtild1**2 = Xprim1 + ! a1 + 2. * a2 * xtild2 + 3. * a3 * xtild2**2 = Xprim2 + + ! (passe par les points (Xf(it), xtild(it)) et (Xf(it + 1), + ! xtild(it + 1)) + + ! On en revient à resoudre un système de 4 équations à 4 inconnues + ! a0, a1, a2, a3. + + use nrtype, only: k8 + + REAL(K8), intent(in):: xf1, xf2, xprim1, xprim2, xtild1, xtild2 + REAL(K8), intent(out):: a0, a1, a2, a3 + + ! Local: + REAL(K8) xtil1car, xtil2car, derr, x1x2car + + !------------------------------------------------------------ + + xtil1car = xtild1 * xtild1 + xtil2car = xtild2 * xtild2 + + derr = 2. * (xf2-xf1)/(xtild1-xtild2) + + x1x2car = (xtild1-xtild2) * (xtild1-xtild2) + + a3 = (derr+xprim1+xprim2)/x1x2car + a2 = (xprim1-xprim2+3. * a3 * (xtil2car-xtil1car))/(2. * (xtild1-xtild2)) + + a1 = xprim1 - 3. * a3 * xtil1car - 2. * a2 * xtild1 + a0 = xf1 - a3 * xtild1 * xtil1car - a2 * xtil1car - a1 * xtild1 + + END SUBROUTINE coefpoly + +end module coefpoly_m Index: LMDZ5/branches/testing/libf/dyn3d_common/defrun.F =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/defrun.F (revision 2220) +++ (revision ) @@ -1,496 +1,0 @@ -! -! $Id$ -! -c -c - SUBROUTINE defrun( tapedef, etatinit, clesphy0 ) -c - USE control_mod - - IMPLICIT NONE -c----------------------------------------------------------------------- -c Auteurs : L. Fairhead , P. Le Van . -c -c Arguments : -c -c tapedef : -c etatinit : = TRUE , on ne compare pas les valeurs des para- -c -metres du zoom avec celles lues sur le fichier start . -c clesphy0 : sortie . -c - LOGICAL etatinit - INTEGER tapedef - - INTEGER longcles - PARAMETER( longcles = 20 ) - REAL clesphy0( longcles ) -c -c Declarations : -c -------------- -#include "dimensions.h" -#include "paramet.h" -#include "logic.h" -#include "serre.h" -#include "comdissnew.h" -#include "clesph0.h" -c -c -c local: -c ------ - - CHARACTER ch1*72,ch2*72,ch3*72,ch4*12 - INTEGER tapeout - REAL clonn,clatt,grossismxx,grossismyy - REAL dzoomxx,dzoomyy,tauxx,tauyy - LOGICAL fxyhypbb, ysinuss - INTEGER i - -c -c ------------------------------------------------------------------- -c -c ......... Version du 29/04/97 .......... -c -c Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot, -c tetatemp ajoutes pour la dissipation . -c -c Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** -c -c Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb. -c Sinon , choix de fxynew , a derivee sinusoidale .. -c -c ...... etatinit = . TRUE. si defrun est appele dans ETAT0_LMD ou -c LIMIT_LMD pour l'initialisation de start.dat (dic) et -c de limit.dat ( dic) ........... -c Sinon etatinit = . FALSE . -c -c Donc etatinit = .F. si on veut comparer les valeurs de grossismx , -c grossismy,clon,clat, fxyhypb lues sur le fichier start avec -c celles passees par run.def , au debut du gcm, apres l'appel a -c lectba . -c Ces parmetres definissant entre autres la grille et doivent etre -c pareils et coherents , sinon il y aura divergence du gcm . -c -c----------------------------------------------------------------------- -c initialisations: -c ---------------- - - tapeout = 6 - -c----------------------------------------------------------------------- -c Parametres de controle du run: -c----------------------------------------------------------------------- - - OPEN( tapedef,file ='gcm.def',status='old',form='formatted') - - - READ (tapedef,9000) ch1,ch2,ch3 - WRITE(tapeout,9000) ch1,ch2,ch3 - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) dayref - WRITE(tapeout,9001) ch1,'dayref' - WRITE(tapeout,*) dayref - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) anneeref - WRITE(tapeout,9001) ch1,'anneeref' - WRITE(tapeout,*) anneeref - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) nday - WRITE(tapeout,9001) ch1,'nday' - WRITE(tapeout,*) nday - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) day_step - WRITE(tapeout,9001) ch1,'day_step' - WRITE(tapeout,*) day_step - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) iperiod - WRITE(tapeout,9001) ch1,'iperiod' - WRITE(tapeout,*) iperiod - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) iapp_tracvl - WRITE(tapeout,9001) ch1,'iapp_tracvl' - WRITE(tapeout,*) iapp_tracvl - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) iconser - WRITE(tapeout,9001) ch1,'iconser' - WRITE(tapeout,*) iconser - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) iecri - WRITE(tapeout,9001) ch1,'iecri' - WRITE(tapeout,*) iecri - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) periodav - WRITE(tapeout,9001) ch1,'periodav' - WRITE(tapeout,*) periodav - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) dissip_period - WRITE(tapeout,9001) ch1,'dissip_period' - WRITE(tapeout,*) dissip_period - -ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... -ccc - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) lstardis - WRITE(tapeout,9001) ch1,'lstardis' - WRITE(tapeout,*) lstardis - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) nitergdiv - WRITE(tapeout,9001) ch1,'nitergdiv' - WRITE(tapeout,*) nitergdiv - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) nitergrot - WRITE(tapeout,9001) ch1,'nitergrot' - WRITE(tapeout,*) nitergrot - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) niterh - WRITE(tapeout,9001) ch1,'niterh' - WRITE(tapeout,*) niterh - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) tetagdiv - WRITE(tapeout,9001) ch1,'tetagdiv' - WRITE(tapeout,*) tetagdiv - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) tetagrot - WRITE(tapeout,9001) ch1,'tetagrot' - WRITE(tapeout,*) tetagrot - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) tetatemp - WRITE(tapeout,9001) ch1,'tetatemp' - WRITE(tapeout,*) tetatemp - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) coefdis - WRITE(tapeout,9001) ch1,'coefdis' - WRITE(tapeout,*) coefdis -c - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) purmats - WRITE(tapeout,9001) ch1,'purmats' - WRITE(tapeout,*) purmats - -c ............................................................... - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) iflag_phys - WRITE(tapeout,9001) ch1,'iflag_phys' - WRITE(tapeout,*) iflag_phys - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) iphysiq - WRITE(tapeout,9001) ch1,'iphysiq' - WRITE(tapeout,*) iphysiq - - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) cycle_diurne - WRITE(tapeout,9001) ch1,'cycle_diurne' - WRITE(tapeout,*) cycle_diurne - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) soil_model - WRITE(tapeout,9001) ch1,'soil_model' - WRITE(tapeout,*) soil_model - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) new_oliq - WRITE(tapeout,9001) ch1,'new_oliq' - WRITE(tapeout,*) new_oliq - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) ok_orodr - WRITE(tapeout,9001) ch1,'ok_orodr' - WRITE(tapeout,*) ok_orodr - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) ok_orolf - WRITE(tapeout,9001) ch1,'ok_orolf' - WRITE(tapeout,*) ok_orolf - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) ok_limitvrai - WRITE(tapeout,9001) ch1,'ok_limitvrai' - WRITE(tapeout,*) ok_limitvrai - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) nbapp_rad - WRITE(tapeout,9001) ch1,'nbapp_rad' - WRITE(tapeout,*) nbapp_rad - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) iflag_con - WRITE(tapeout,9001) ch1,'iflag_con' - WRITE(tapeout,*) iflag_con - - DO i = 1, longcles - clesphy0(i) = 0. - ENDDO - clesphy0(1) = REAL( iflag_con ) - clesphy0(2) = REAL( nbapp_rad ) - - IF( cycle_diurne ) clesphy0(3) = 1. - IF( soil_model ) clesphy0(4) = 1. - IF( new_oliq ) clesphy0(5) = 1. - IF( ok_orodr ) clesphy0(6) = 1. - IF( ok_orolf ) clesphy0(7) = 1. - IF( ok_limitvrai ) clesphy0(8) = 1. - - -ccc .... P. Le Van , ajout le 7/03/95 .pour le zoom ... -c ......... ( modif le 17/04/96 ) ......... -c - IF( etatinit ) GO TO 100 - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) clonn - WRITE(tapeout,9001) ch1,'clon' - WRITE(tapeout,*) clonn - IF( ABS(clon - clonn).GE. 0.001 ) THEN - WRITE(tapeout,*) ' La valeur de clon passee par run.def est diffe - *rente de celle lue sur le fichier start ' - STOP - ENDIF -c - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) clatt - WRITE(tapeout,9001) ch1,'clat' - WRITE(tapeout,*) clatt - - IF( ABS(clat - clatt).GE. 0.001 ) THEN - WRITE(tapeout,*) ' La valeur de clat passee par run.def est diffe - *rente de celle lue sur le fichier start ' - STOP - ENDIF - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) grossismxx - WRITE(tapeout,9001) ch1,'grossismx' - WRITE(tapeout,*) grossismxx - - IF( ABS(grossismx - grossismxx).GE. 0.001 ) THEN - WRITE(tapeout,*) ' La valeur de grossismx passee par run.def est - , differente de celle lue sur le fichier start ' - STOP - ENDIF - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) grossismyy - WRITE(tapeout,9001) ch1,'grossismy' - WRITE(tapeout,*) grossismyy - - IF( ABS(grossismy - grossismyy).GE. 0.001 ) THEN - WRITE(tapeout,*) ' La valeur de grossismy passee par run.def est - , differente de celle lue sur le fichier start ' - STOP - ENDIF - - IF( grossismx.LT.1. ) THEN - WRITE(tapeout,*) ' *** ATTENTION !! grossismx < 1 . *** ' - STOP - ELSE - alphax = 1. - 1./ grossismx - ENDIF - - - IF( grossismy.LT.1. ) THEN - WRITE(tapeout,*) ' *** ATTENTION !! grossismy < 1 . *** ' - STOP - ELSE - alphay = 1. - 1./ grossismy - ENDIF - -c -c alphax et alphay sont les anciennes formulat. des grossissements -c -c - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) fxyhypbb - WRITE(tapeout,9001) ch1,'fxyhypbb' - WRITE(tapeout,*) fxyhypbb - - IF( .NOT.fxyhypb ) THEN - IF( fxyhypbb ) THEN - WRITE(tapeout,*) ' ******** PBS DANS DEFRUN ******** ' - WRITE(tapeout,*)' *** fxyhypb lu sur le fichier start est F' - *, ' alors qu il est T sur run.def ***' - STOP - ENDIF - ELSE - IF( .NOT.fxyhypbb ) THEN - WRITE(tapeout,*) ' ******** PBS DANS DEFRUN ******** ' - WRITE(tapeout,*)' *** fxyhypb lu sur le fichier start est t' - *, ' alors qu il est F sur run.def ***' - STOP - ENDIF - ENDIF -c - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) dzoomxx - WRITE(tapeout,9001) ch1,'dzoomx' - WRITE(tapeout,*) dzoomxx - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) dzoomyy - WRITE(tapeout,9001) ch1,'dzoomy' - WRITE(tapeout,*) dzoomyy - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) tauxx - WRITE(tapeout,9001) ch1,'taux' - WRITE(tapeout,*) tauxx - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) tauyy - WRITE(tapeout,9001) ch1,'tauy' - WRITE(tapeout,*) tauyy - - IF( fxyhypb ) THEN - - IF( ABS(dzoomx - dzoomxx).GE. 0.001 ) THEN - WRITE(tapeout,*)' La valeur de dzoomx passee par run.def est dif - *ferente de celle lue sur le fichier start ' - CALL ABORT_gcm("defrun", "", 1) - ENDIF - - IF( ABS(dzoomy - dzoomyy).GE. 0.001 ) THEN - WRITE(tapeout,*)' La valeur de dzoomy passee par run.def est dif - *ferente de celle lue sur le fichier start ' - CALL ABORT_gcm("defrun", "", 1) - ENDIF - - IF( ABS(taux - tauxx).GE. 0.001 ) THEN - WRITE(6,*)' La valeur de taux passee par run.def est differente - * de celle lue sur le fichier start ' - CALL ABORT_gcm("defrun", "", 1) - ENDIF - - IF( ABS(tauy - tauyy).GE. 0.001 ) THEN - WRITE(6,*)' La valeur de tauy passee par run.def est differente - * de celle lue sur le fichier start ' - CALL ABORT_gcm("defrun", "", 1) - ENDIF - - ENDIF - -cc - IF( .NOT.fxyhypb ) THEN - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) ysinuss - WRITE(tapeout,9001) ch1,'ysinus' - WRITE(tapeout,*) ysinuss - - - IF( .NOT.ysinus ) THEN - IF( ysinuss ) THEN - WRITE(6,*) ' ******** PBS DANS DEFRUN ******** ' - WRITE(tapeout,*)'** ysinus lu sur le fichier start est F', - * ' alors qu il est T sur run.def ***' - STOP - ENDIF - ELSE - IF( .NOT.ysinuss ) THEN - WRITE(6,*) ' ******** PBS DANS DEFRUN ******** ' - WRITE(tapeout,*)'** ysinus lu sur le fichier start est T', - * ' alors qu il est F sur run.def ***' - STOP - ENDIF - ENDIF - ENDIF -c - WRITE(6,*) ' alphax alphay defrun ',alphax,alphay - - CLOSE(tapedef) - - RETURN -c ............................................... -c -100 CONTINUE -c - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) clon - WRITE(tapeout,9001) ch1,'clon' - WRITE(tapeout,*) clon -c - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) clat - WRITE(tapeout,9001) ch1,'clat' - WRITE(tapeout,*) clat - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) grossismx - WRITE(tapeout,9001) ch1,'grossismx' - WRITE(tapeout,*) grossismx - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) grossismy - WRITE(tapeout,9001) ch1,'grossismy' - WRITE(tapeout,*) grossismy - - IF( grossismx.LT.1. ) THEN - WRITE(tapeout,*) '*** ATTENTION !! grossismx < 1 . *** ' - STOP - ELSE - alphax = 1. - 1./ grossismx - ENDIF - - IF( grossismy.LT.1. ) THEN - WRITE(tapeout,*) ' *** ATTENTION !! grossismy < 1 . *** ' - STOP - ELSE - alphay = 1. - 1./ grossismy - ENDIF - -c - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) fxyhypb - WRITE(tapeout,9001) ch1,'fxyhypb' - WRITE(tapeout,*) fxyhypb - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) dzoomx - WRITE(tapeout,9001) ch1,'dzoomx' - WRITE(tapeout,*) dzoomx - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) dzoomy - WRITE(tapeout,9001) ch1,'dzoomy' - WRITE(tapeout,*) dzoomy - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) taux - WRITE(tapeout,9001) ch1,'taux' - WRITE(tapeout,*) taux -c - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) tauy - WRITE(tapeout,9001) ch1,'tauy' - WRITE(tapeout,*) tauy - - READ (tapedef,9001) ch1,ch4 - READ (tapedef,*) ysinus - WRITE(tapeout,9001) ch1,'ysinus' - WRITE(tapeout,*) ysinus - - WRITE(tapeout,*) ' alphax alphay defrun ',alphax,alphay -c -9000 FORMAT(3(/,a72)) -9001 FORMAT(/,a72,/,a12) -cc - CLOSE(tapedef) - - RETURN - END Index: LMDZ5/branches/testing/libf/dyn3d_common/fxhyp_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/fxhyp_m.F90 (revision 2258) +++ LMDZ5/branches/testing/libf/dyn3d_common/fxhyp_m.F90 (revision 2258) @@ -0,0 +1,250 @@ +module fxhyp_m + + IMPLICIT NONE + +contains + + SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) + + ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32 + ! Author: P. Le Van, from formulas by R. Sadourny + + ! Calcule les longitudes et dérivées dans la grille du GCM pour + ! une fonction f(x) à dérivée tangente hyperbolique. + + ! Il vaut mieux avoir : grossismx \times dzoom < pi + + ! Le premier point scalaire pour une grille regulière (grossismx = + ! 1., taux=0., clon=0.) est à - 180 degrés. + + use arth_m, only: arth + use invert_zoom_x_m, only: invert_zoom_x, nmax + use nrtype, only: pi, pi_d, twopi, twopi_d, k8 + use principal_cshift_m, only: principal_cshift + + include "dimensions.h" + ! for iim + + include "serre.h" + ! for clon, grossismx, dzoomx, taux + + REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1) + real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1) + + ! Local: + real rlonm025(iim + 1), rlonp025(iim + 1) + REAL dzoom, step + real d_rlonv(iim) + REAL(K8) xtild(0:2 * nmax) + REAL(K8) fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax) + REAL(K8) Xf(0:2 * nmax), xxpr(2 * nmax) + REAL(K8) fa, fb + INTEGER i, is2 + REAL(K8) xmoy, fxm + + !---------------------------------------------------------------------- + + print *, "Call sequence information: fxhyp" + + test_iim: if (iim==1) then + rlonv(1)=0. + rlonu(1)=pi + rlonv(2)=rlonv(1)+twopi + rlonu(2)=rlonu(1)+twopi + + xprimm025(:)=1. + xprimv(:)=1. + xprimu(:)=1. + xprimp025(:)=1. + else test_iim + test_grossismx: if (grossismx == 1.) then + step = twopi / iim + + xprimm025(:iim) = step + xprimp025(:iim) = step + xprimv(:iim) = step + xprimu(:iim) = step + + rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim) + rlonm025(:iim) = rlonv(:iim) - 0.25 * step + rlonp025(:iim) = rlonv(:iim) + 0.25 * step + rlonu(:iim) = rlonv(:iim) + 0.5 * step + else test_grossismx + dzoom = dzoomx * twopi_d + xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1) + + ! Compute fhyp: + DO i = nmax, 2 * nmax + fa = taux * (dzoom / 2. - xtild(i)) + fb = xtild(i) * (pi_d - xtild(i)) + + IF (200. * fb < - fa) THEN + fhyp(i) = - 1. + ELSE IF (200. * fb < fa) THEN + fhyp(i) = 1. + ELSE + IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN + IF (200. * fb + fa < 1e-10) THEN + fhyp(i) = - 1. + ELSE IF (200. * fb - fa < 1e-10) THEN + fhyp(i) = 1. + END IF + ELSE + fhyp(i) = TANH(fa / fb) + END IF + END IF + + IF (xtild(i) == 0.) fhyp(i) = 1. + IF (xtild(i) == pi_d) fhyp(i) = -1. + END DO + + ! Calcul de beta + + ffdx = 0. + + DO i = nmax + 1, 2 * nmax + xmoy = 0.5 * (xtild(i-1) + xtild(i)) + fa = taux * (dzoom / 2. - xmoy) + fb = xmoy * (pi_d - xmoy) + + IF (200. * fb < - fa) THEN + fxm = - 1. + ELSE IF (200. * fb < fa) THEN + fxm = 1. + ELSE + IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN + IF (200. * fb + fa < 1e-10) THEN + fxm = - 1. + ELSE IF (200. * fb - fa < 1e-10) THEN + fxm = 1. + END IF + ELSE + fxm = TANH(fa / fb) + END IF + END IF + + IF (xmoy == 0.) fxm = 1. + IF (xmoy == pi_d) fxm = -1. + + ffdx = ffdx + fxm * (xtild(i) - xtild(i-1)) + END DO + + print *, "ffdx = ", ffdx + beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d) + print *, "beta = ", beta + + IF (2. * beta - grossismx <= 0.) THEN + print *, 'Bad choice of grossismx, taux, dzoomx.' + print *, 'Decrease dzoomx or grossismx.' + STOP 1 + END IF + + ! calcul de Xprimt + Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp + xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1) + + ! Calcul de Xf + + DO i = nmax + 1, 2 * nmax + xmoy = 0.5 * (xtild(i-1) + xtild(i)) + fa = taux * (dzoom / 2. - xmoy) + fb = xmoy * (pi_d - xmoy) + + IF (200. * fb < - fa) THEN + fxm = - 1. + ELSE IF (200. * fb < fa) THEN + fxm = 1. + ELSE + fxm = TANH(fa / fb) + END IF + + IF (xmoy == 0.) fxm = 1. + IF (xmoy == pi_d) fxm = -1. + xxpr(i) = beta + (grossismx - beta) * fxm + END DO + + xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1) + + Xf(0) = - pi_d + + DO i=1, 2 * nmax - 1 + Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1)) + END DO + + Xf(2 * nmax) = pi_d + + call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), & + xprimm025(:iim), xuv = - 0.25_k8) + call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), & + xuv = 0._k8) + call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), & + xuv = 0.5_k8) + call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), & + xprimp025(:iim), xuv = 0.25_k8) + end if test_grossismx + + 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 + + if (rlonm025(is2) < - pi) then + print *, 'Rlonm025 plus petit que - pi !' + STOP 1 + end if + 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) * 180. / pi, & + "degrees" + print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, & + "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 if test_iim + + END SUBROUTINE fxhyp + +end module fxhyp_m Index: LMDZ5/branches/testing/libf/dyn3d_common/fxyhyper.F =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/fxyhyper.F (revision 2220) +++ (revision ) @@ -1,139 +1,0 @@ -! -! $Header$ -! -c -c - SUBROUTINE fxyhyper ( yzoom, grossy, dzoomy,tauy , - , xzoom, grossx, dzoomx,taux , - , rlatu,yprimu,rlatv,yprimv,rlatu1, yprimu1, rlatu2, yprimu2 , - , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025) - - IMPLICIT NONE -c -c Auteur : P. Le Van . -c -c d'apres formulations de R. Sadourny . -c -c -c Ce spg calcule les latitudes( routine fyhyp ) et longitudes( fxhyp ) -c par des fonctions a tangente hyperbolique . -c -c Il y a 3 parametres ,en plus des coordonnees du centre du zoom (xzoom -c et yzoom ) : -c -c a) le grossissement du zoom : grossy ( en y ) et grossx ( en x ) -c b) l' extension du zoom : dzoomy ( en y ) et dzoomx ( en x ) -c c) la raideur de la transition du zoom : taux et tauy -c -c N.B : Il vaut mieux avoir : grossx * dzoomx < pi ( radians ) -c ****** -c et grossy * dzoomy < pi/2 ( radians ) -c -#include "dimensions.h" -#include "paramet.h" - - -c ..... Arguments ... -c - REAL xzoom,yzoom,grossx,grossy,dzoomx,dzoomy,taux,tauy - REAL rlatu(jjp1), yprimu(jjp1),rlatv(jjm), yprimv(jjm), - , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm) - REAL rlonu(iip1),xprimu(iip1),rlonv(iip1),xprimv(iip1), - , rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),xprimp025(iip1) - REAL(KIND=8) dxmin, dxmax , dymin, dymax - -c .... var. locales ..... -c - INTEGER i,j -c - - CALL fyhyp ( yzoom, grossy, dzoomy,tauy , - , rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 , - , dymin,dymax ) - - CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv, - , xprimv,rlonu,xprimu,rlonp025,xprimp025 , dxmin,dxmax ) - - - DO i = 1, iip1 - IF(rlonp025(i).LT.rlonv(i)) THEN - WRITE(6,*) ' Attention ! rlonp025 < rlonv',i - STOP - ENDIF - - IF(rlonv(i).LT.rlonm025(i)) THEN - WRITE(6,*) ' Attention ! rlonm025 > rlonv',i - STOP - ENDIF - - IF(rlonp025(i).GT.rlonu(i)) THEN - WRITE(6,*) ' Attention ! rlonp025 > rlonu',i - STOP - ENDIF - ENDDO - - WRITE(6,*) ' *** TEST DE COHERENCE OK POUR FX **** ' - -c - DO j = 1, jjm -c - IF(rlatu1(j).LE.rlatu2(j)) THEN - WRITE(6,*)'Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j - STOP 13 - ENDIF -c - IF(rlatu2(j).LE.rlatu(j+1)) THEN - WRITE(6,*)'Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j - STOP 14 - ENDIF -c - IF(rlatu(j).LE.rlatu1(j)) THEN - WRITE(6,*)' Attention ! rlatu < rlatu1 ',rlatu(j),rlatu1(j),j - STOP 15 - ENDIF -c - IF(rlatv(j).LE.rlatu2(j)) THEN - WRITE(6,*)' Attention ! rlatv < rlatu2 ',rlatv(j),rlatu2(j),j - STOP 16 - ENDIF -c - IF(rlatv(j).ge.rlatu1(j)) THEN - WRITE(6,*)' Attention ! rlatv > rlatu1 ',rlatv(j),rlatu1(j),j - STOP 17 - ENDIF -c - IF(rlatv(j).ge.rlatu(j)) THEN - WRITE(6,*) ' Attention ! rlatv > rlatu ',rlatv(j),rlatu(j),j - STOP 18 - ENDIF -c - ENDDO -c - WRITE(6,*) ' *** TEST DE COHERENCE OK POUR FY **** ' -c - WRITE(6,18) - WRITE(6,*) ' Latitudes ' - WRITE(6,*) ' *********** ' - WRITE(6,18) - WRITE(6,3) dymin, dymax - WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par - ,ametres grossism , tau , dzoom pour Y et repasser ! ' -c - WRITE(6,18) - WRITE(6,*) ' Longitudes ' - WRITE(6,*) ' ************ ' - WRITE(6,18) - WRITE(6,3) dxmin, dxmax - WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par - ,ametres grossism , tau , dzoom pour Y et repasser ! ' - WRITE(6,18) -c -3 Format(1x, ' Au centre du zoom , la longueur de la maille est', - , ' d environ ',f8.2 ,' degres ', - , ' alors que la maille en dehors de la zone du zoom est d environ - , ', f8.2,' degres ' ) -18 FORMAT(/) - - RETURN - END - Index: LMDZ5/branches/testing/libf/dyn3d_common/fyhyp.F =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/fyhyp.F (revision 2220) +++ (revision ) @@ -1,378 +1,0 @@ -! -! $Id$ -! -c -c - SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau , - , rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 , - , champmin,champmax ) - -cc ... Version du 01/04/2001 .... - - IMPLICIT NONE -c -c ... Auteur : P. Le Van ... -c -c ....... d'apres formulations de R. Sadourny ....... -c -c Calcule les latitudes et derivees dans la grille du GCM pour une -c fonction f(y) a tangente hyperbolique . -c -c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc) -c dzoom etant la distance totale de la zone du zoom ( en radians ) -c tau la raideur de la transition de l'interieur a l'exterieur du zoom -c -c -c N.B : Il vaut mieux avoir : grossism * dzoom < pi/2 (radians) ,en lati. -c ******************************************************************** -c -c -#include "dimensions.h" -#include "paramet.h" - - INTEGER nmax , nmax2 - PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) -c -c -c ....... arguments d'entree ....... -c - REAL yzoomdeg, grossism,dzooma,tau -c ( rentres par run.def ) - -c ....... arguments de sortie ....... -c - REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm), - , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm) - -c -c ..... champs locaux ..... -c - - REAL dzoom - REAL(KIND=8) ylat(jjp1), yprim(jjp1) - REAL(KIND=8) yuv - REAL(KIND=8) yt(0:nmax2) - REAL(KIND=8) fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2) - SAVE Ytprim, yt,Yf - REAL(KIND=8) Yf(0:nmax2),yypr(0:nmax2) - REAL(KIND=8) yvrai(jjp1), yprimm(jjp1),ylatt(jjp1) - REAL(KIND=8) pi,depi,pis2,epsilon,y0,pisjm - REAL(KIND=8) yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax - REAL(KIND=8) yfi,Yf1,ffdy - REAL(KIND=8) ypn,deply,y00 - SAVE y00, deply - - INTEGER i,j,it,ik,iter,jlat - INTEGER jpn,jjpn - SAVE jpn - REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m - REAL(KIND=8) fa(0:nmax2),fb(0:nmax2) - REAL y0min,y0max - - REAL(KIND=8) heavyside - - pi = 2. * ASIN(1.) - depi = 2. * pi - pis2 = pi/2. - pisjm = pi/ REAL(jjm) - epsilon = 1.e-3 - y0 = yzoomdeg * pi/180. - - IF( dzooma.LT.1.) THEN - dzoom = dzooma * pi - ELSEIF( dzooma.LT. 12. ) THEN - WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug - ,menter et relancer ! ' - STOP 1 - ELSE - dzoom = dzooma * pi/180. - ENDIF - - WRITE(6,18) - WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)' - WRITE(6,24) y0,grossism,tau,dzoom - - DO i = 0, nmax2 - yt(i) = - pis2 + REAL(i)* pi /nmax2 - ENDDO - - 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).LT.y0 ) THEN - fa (i) = tau* (yt(i)-y0+dzoom/2. ) - fb(i) = (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) ) - ELSEIF ( yt(i).GT.y0 ) THEN - fa(i) = tau *(y0-yt(i)+dzoom/2. ) - fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 ) - ENDIF - - IF( 200.* fb(i) .LT. - fa(i) ) THEN - fhyp ( i) = - 1. - ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN - fhyp ( i) = 1. - ELSE - fhyp(i) = TANH ( fa(i)/fb(i) ) - ENDIF - - IF( yt(i).EQ.y0 ) fhyp(i) = 1. - IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1. - - ENDDO - -cc .... Calcul de beta .... -c - ffdy = 0. - - DO i = 1, nmax2 - ymoy = 0.5 * ( yt(i-1) + yt( i ) ) - IF( ymoy.LT.y0 ) THEN - fa(i)= tau * ( ymoy-y0+dzoom/2.) - fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy ) - ELSEIF ( ymoy.GT.y0 ) THEN - fa(i)= tau * ( y0-ymoy+dzoom/2. ) - fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 ) - ENDIF - - IF( 200.* fb(i) .LT. - fa(i) ) THEN - fxm ( i) = - 1. - ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN - fxm ( i) = 1. - ELSE - fxm(i) = TANH ( fa(i)/fb(i) ) - ENDIF - IF( ymoy.EQ.y0 ) fxm(i) = 1. - IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1. - ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) ) - - ENDDO - - beta = ( grossism * ffdy - pi ) / ( ffdy - pi ) - - IF( 2.*beta - grossism.LE. 0.) THEN - - WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou - ,tine fyhyp est mauvaise ! ' - WRITE(6,*)'Modifier les valeurs de grossismy ,tauy ou dzoomy', - , ' et relancer ! *** ' - CALL ABORT_GCM("FYHYP", "", 1) - - ENDIF -c -c ..... calcul de Ytprim ..... -c - - DO i = 0, nmax2 - Ytprim(i) = beta + ( grossism - beta ) * fhyp(i) - ENDDO - -c ..... Calcul de Yf ........ - - Yf(0) = - pis2 - DO i = 1, nmax2 - yypr(i) = beta + ( grossism - beta ) * fxm(i) - ENDDO - - DO i=1,nmax2 - Yf(i) = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) ) - ENDDO - -c **************************************************************** -c -c ..... yuv = 0. si calcul des latitudes aux pts. U ..... -c ..... yuv = 0.5 si calcul des latitudes aux pts. V ..... -c - WRITE(6,18) -c - DO 5000 ik = 1,4 - - IF( ik.EQ.1 ) THEN - yuv = 0. - jlat = jjm + 1 - ELSE IF ( ik.EQ.2 ) THEN - yuv = 0.5 - jlat = jjm - ELSE IF ( ik.EQ.3 ) THEN - yuv = 0.25 - jlat = jjm - ELSE IF ( ik.EQ.4 ) THEN - yuv = 0.75 - jlat = jjm - ENDIF -c - yo1 = 0. - DO 1500 j = 1,jlat - yo1 = 0. - ylon2 = - pis2 + pisjm * ( REAL(j) + yuv -1.) - yfi = ylon2 -c - DO 250 it = nmax2,0,-1 - IF( yfi.GE.Yf(it)) GO TO 350 -250 CONTINUE - it = 0 -350 CONTINUE - - yi = yt(it) - IF(it.EQ.nmax2) THEN - it = nmax2 -1 - Yf(it+1) = pis2 - ENDIF -c ................................................................. -c .... Interpolation entre yi(it) et yi(it+1) pour avoir Y(yi) -c ..... et Y'(yi) ..... -c ................................................................. - - 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 - - DO 500 iter = 1,300 - yi = yi - ( Yf1 - yfi )/ Yprimin - - IF( ABS(yi-yo1).LE.epsilon) GO TO 550 - yo1 = yi - yi2 = yi * yi - Yf1 = a0 + a1 * yi + a2 * yi2 + a3 * yi2 * yi - Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi2 -500 CONTINUE - WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter - STOP 2 -550 CONTINUE -c - Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi* yi - yprim(j) = pi / ( jjm * Yprimin ) - yvrai(j) = yi - -1500 CONTINUE - - DO j = 1, jlat -1 - IF( yvrai(j+1). LT. yvrai(j) ) THEN - WRITE(6,*) ' PBS. avec rlat(',j+1,') plus petit que rlat(',j, - , ')' - STOP 3 - ENDIF - ENDDO - - WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2' - , ,' et pi/2 ' -c - IF( ik.EQ.1 ) THEN - ypn = pis2 - DO j = jlat,1,-1 - IF( yvrai(j).LE. ypn ) GO TO 1502 - ENDDO -1502 CONTINUE - - jpn = j - y00 = yvrai(jpn) - deply = pis2 - y00 - ENDIF - - DO j = 1, jjm +1 - jpn - ylatt (j) = -pis2 - y00 + yvrai(jpn+j-1) - yprimm(j) = yprim(jpn+j-1) - ENDDO - - jjpn = jpn - IF( jlat.EQ. jjm ) jjpn = jpn -1 - - DO j = 1,jjpn - ylatt (j + jjm+1 -jpn) = yvrai(j) + deply - yprimm(j + jjm+1 -jpn) = yprim(j) - ENDDO - -c *********** Fin de la reorganisation ************* -c - 1600 CONTINUE - - DO j = 1, jlat - ylat(j) = ylatt( jlat +1 -j ) - yprim(j) = yprimm( jlat +1 -j ) - ENDDO - - DO j = 1, jlat - yvrai(j) = ylat(j)*180./pi - ENDDO - - IF( ik.EQ.1 ) THEN -c WRITE(6,18) -c WRITE(6,*) ' YLAT en U apres ( en deg. ) ' -c WRITE(6,68) (yvrai(j),j=1,jlat) -cc WRITE(6,*) ' YPRIM ' -cc WRITE(6,445) ( yprim(j),j=1,jlat) - - DO j = 1, jlat - rrlatu(j) = ylat( j ) - yyprimu(j) = yprim( j ) - ENDDO - - ELSE IF ( ik.EQ. 2 ) THEN -c WRITE(6,18) -c WRITE(6,*) ' YLAT en V apres ( en deg. ) ' -c WRITE(6,68) (yvrai(j),j=1,jlat) -cc WRITE(6,*)' YPRIM ' -cc WRITE(6,445) ( yprim(j),j=1,jlat) - - DO j = 1, jlat - rrlatv(j) = ylat( j ) - yyprimv(j) = yprim( j ) - ENDDO - - ELSE IF ( ik.EQ. 3 ) THEN -c WRITE(6,18) -c WRITE(6,*) ' YLAT en U + 0.75 apres ( en deg. ) ' -c WRITE(6,68) (yvrai(j),j=1,jlat) -cc WRITE(6,*) ' YPRIM ' -cc WRITE(6,445) ( yprim(j),j=1,jlat) - - DO j = 1, jlat - rlatu2(j) = ylat( j ) - yprimu2(j) = yprim( j ) - ENDDO - - ELSE IF ( ik.EQ. 4 ) THEN -c WRITE(6,18) -c WRITE(6,*) ' YLAT en U + 0.25 apres ( en deg. ) ' -c WRITE(6,68)(yvrai(j),j=1,jlat) -cc WRITE(6,*) ' YPRIM ' -cc WRITE(6,68) ( yprim(j),j=1,jlat) - - DO j = 1, jlat - rlatu1(j) = ylat( j ) - yprimu1(j) = yprim( j ) - ENDDO - - ENDIF - -5000 CONTINUE -c - WRITE(6,18) -c -c ..... fin de la boucle do 5000 ..... - - DO j = 1, jjm - ylat(j) = rrlatu(j) - rrlatu(j+1) - ENDDO - champmin = 1.e12 - champmax = -1.e12 - DO j = 1, jjm - champmin = MIN( champmin, ylat(j) ) - champmax = MAX( champmax, ylat(j) ) - ENDDO - champmin = champmin * 180./pi - champmax = champmax * 180./pi - -24 FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3) -18 FORMAT(/) -68 FORMAT(1x,7f9.2) - - RETURN - END Index: LMDZ5/branches/testing/libf/dyn3d_common/fyhyp_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/fyhyp_m.F90 (revision 2258) +++ LMDZ5/branches/testing/libf/dyn3d_common/fyhyp_m.F90 (revision 2258) @@ -0,0 +1,342 @@ +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 + use nrtype, only: k8 + + 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: + + REAL(K8) champmin, champmax + INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax + REAL dzoom ! distance totale de la zone du zoom (en radians) + REAL(K8) ylat(jjm + 1), yprim(jjm + 1) + REAL(K8) yuv + REAL(K8), save:: yt(0:nmax2) + REAL(K8) fhyp(0:nmax2), beta + REAL(K8), save:: ytprim(0:nmax2) + REAL(K8) fxm(0:nmax2) + REAL(K8), save:: yf(0:nmax2) + REAL(K8) yypr(0:nmax2) + REAL(K8) yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1) + REAL(K8) pi, pis2, epsilon, y0, pisjm + REAL(K8) yo1, yi, ylon2, ymoy, yprimin + REAL(K8) yfi, yf1, ffdy + REAL(K8) ypn, deply, y00 + SAVE y00, deply + + INTEGER i, j, it, ik, iter, jlat + INTEGER jpn, jjpn + SAVE jpn + REAL(K8) a0, a1, a2, a3, yi2, heavyy0, heavyy0m + REAL(K8) fa(0:nmax2), fb(0:nmax2) + REAL y0min, y0max + + REAL(K8) 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 Index: LMDZ5/branches/testing/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90 =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90 (revision 2220) +++ LMDZ5/branches/testing/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90 (revision 2258) @@ -31,7 +31,4 @@ INTEGER out_latudim,out_latvdim,out_dim(3) INTEGER out_levdim - - INTEGER, PARAMETER :: longcles = 20 - REAL clesphy0(longcles) INTEGER start(4),COUNT(4) @@ -60,5 +57,5 @@ pa= 50000. - CALL conf_gcm( 99, .TRUE. , clesphy0 ) + CALL conf_gcm( 99, .TRUE. ) CALL iniconst CALL inigeom Index: LMDZ5/branches/testing/libf/dyn3d_common/inigeom.F =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/inigeom.F (revision 2220) +++ LMDZ5/branches/testing/libf/dyn3d_common/inigeom.F (revision 2258) @@ -16,4 +16,6 @@ c c + use fxhyp_m, only: fxhyp + use fyhyp_m, only: fyhyp IMPLICIT NONE c @@ -264,9 +266,6 @@ WRITE(6,*)'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' - CALL fxyhyper( clat, grossismy, dzoomy, tauy , - , clon, grossismx, dzoomx, taux , - , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 , - , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 ) - + CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1) + CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) ENDIF Index: LMDZ5/branches/testing/libf/dyn3d_common/invert_zoom_x_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/invert_zoom_x_m.F90 (revision 2258) +++ LMDZ5/branches/testing/libf/dyn3d_common/invert_zoom_x_m.F90 (revision 2258) @@ -0,0 +1,90 @@ +module invert_zoom_x_m + + implicit none + + INTEGER, PARAMETER:: nmax = 30000 + +contains + + subroutine invert_zoom_x(xf, xtild, Xprimt, xlon, xprimm, xuv) + + use coefpoly_m, only: coefpoly + use nrtype, only: pi, pi_d, twopi_d, k8 + + include "dimensions.h" + ! for iim + + include "serre.h" + ! for clon + + REAL(K8), intent(in):: Xf(0:), xtild(0:), Xprimt(0:) ! (0:2 * nmax) + real, intent(out):: xlon(:), xprimm(:) ! (iim) + + REAL(K8), intent(in):: xuv + ! 0. si calcul aux points scalaires + ! 0.5 si calcul aux points U + + ! Local: + REAL(K8) xo1, Xfi, a0, a1, a2, a3, Xf1, Xprimin + integer i, it, iter + REAL(K8), parameter:: my_eps = 1e-6_k8 + + REAL(K8) xxprim(iim), xvrai(iim) + ! intermediary variables because xlon and xprimm are simple precision + + !------------------------------------------------------------------ + + DO i = 1, iim + Xfi = - pi_d + (i + xuv - 0.75_k8) * twopi_d / iim + + it = 2 * nmax + do while (xfi < xf(it) .and. it >= 1) + it = it - 1 + end do + + ! Calcul de Xf(xvrai(i)) + + xvrai(i) = xtild(it) + + IF (it == 2 * nmax) THEN + it = 2 * nmax -1 + END IF + + CALL coefpoly(Xf(it), Xf(it + 1), Xprimt(it), Xprimt(it + 1), & + xtild(it), xtild(it + 1), a0, a1, a2, a3) + Xf1 = Xf(it) + Xprimin = a1 + xvrai(i) * (2._k8 * a2 + xvrai(i) * 3._k8 * a3) + xo1 = xvrai(i) + iter = 1 + + do + xvrai(i) = xvrai(i) - (Xf1 - Xfi) / Xprimin + IF (ABS(xvrai(i) - xo1) <= my_eps .or. iter == 300) exit + xo1 = xvrai(i) + Xf1 = a0 + xvrai(i) * (a1 + xvrai(i) * (a2 + xvrai(i) * a3)) + Xprimin = a1 + xvrai(i) * (2._k8 * a2 + xvrai(i) * 3._k8 * a3) + end DO + + if (ABS(xvrai(i) - xo1) > my_eps) then + ! iter == 300 + print *, 'Pas de solution.' + print *, i, xfi + STOP 1 + end if + + xxprim(i) = twopi_d / (iim * Xprimin) + 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 / 180. * pi + xprimm = xxprim + + end subroutine invert_zoom_x + +end module invert_zoom_x_m Index: LMDZ5/branches/testing/libf/dyn3d_common/principal_cshift_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/dyn3d_common/principal_cshift_m.F90 (revision 2258) +++ LMDZ5/branches/testing/libf/dyn3d_common/principal_cshift_m.F90 (revision 2258) @@ -0,0 +1,45 @@ +module principal_cshift_m + + implicit none + +contains + + 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 nrtype, only: twopi + + include "dimensions.h" + ! for iim + + include "serre.h" + ! for clon + + 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 + +end module principal_cshift_m