Changeset 1441
- Timestamp:
- Jun 4, 2015, 10:21:20 AM (9 years ago)
- Location:
- trunk
- Files:
-
- 3 added
- 15 edited
- 10 moved
-
DOC/chantiers/commit_importants.log (modified) (1 diff)
-
LMDZ.COMMON/libf/dyn3d/advtrac.F90 (modified) (1 diff)
-
LMDZ.COMMON/libf/dyn3d/conf_gcm.F90 (modified) (4 diffs)
-
LMDZ.COMMON/libf/dyn3d/fluxstokenc.F (modified) (1 diff)
-
LMDZ.COMMON/libf/dyn3d/gcm.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d/gcm.F) (5 diffs)
-
LMDZ.COMMON/libf/dyn3d/guide_mod.F90 (modified) (6 diffs)
-
LMDZ.COMMON/libf/dyn3d_common/fxhyp_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/fxhyp.F) (1 diff)
-
LMDZ.COMMON/libf/dyn3d_common/fyhyp_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/fyhyp.F) (1 diff)
-
LMDZ.COMMON/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90 (modified) (2 diffs)
-
LMDZ.COMMON/libf/dyn3d_common/inigeom.F (modified) (2 diffs)
-
LMDZ.COMMON/libf/dyn3d_common/invert_zoom_x_m.F90 (added)
-
LMDZ.COMMON/libf/dyn3d_common/principal_cshift_m.F90 (added)
-
LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90 (modified) (1 diff)
-
LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F90 (modified) (4 diffs)
-
LMDZ.COMMON/libf/dyn3dpar/covcont_p.F (modified) (1 diff)
-
LMDZ.COMMON/libf/dyn3dpar/mod_hallo.F90 (modified) (1 diff)
-
LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90 (modified) (4 diffs)
-
LMDZ.COMMON/libf/grid/dimension/makdim (modified) (1 diff)
-
LMDZ.COMMON/libf/misc/arth_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/misc/arth.F90)
-
LMDZ.COMMON/libf/misc/coefpoly_m.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/coefpoly.F) (1 diff)
-
LMDZ.COMMON/libf/misc/juldate.F (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/juldate.F)
-
LMDZ.COMMON/libf/misc/minmax.F (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/minmax.F)
-
LMDZ.COMMON/libf/misc/minmax2.F (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/minmax2.F)
-
LMDZ.COMMON/libf/misc/nrtype.F90 (added)
-
LMDZ.COMMON/libf/misc/ran1.F (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/ran1.F)
-
LMDZ.COMMON/libf/misc/sort.F (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/sort.F)
-
LMDZ.COMMON/libf/misc/wxios.F90 (modified) (25 diffs)
-
LMDZ.COMMON/makelmdz_fcm (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/DOC/chantiers/commit_importants.log
r1403 r1441 1523 1523 * phy***/dyn1d: this subdirectory contains the 1D model using physics 1524 1524 from phy*** 1525 1526 ********************** 1527 **** commit_v1441 **** 1528 ********************** 1529 Ehouarn: Updates in the dynamics (seq and //) to keep up with updates 1530 in LMDZ5 (up to LMDZ5 trunk, rev 2250): 1531 * compilation: 1532 - added test in grid/dimension/makdim to check that # of longitudes 1533 is a multiple of 8 1534 1535 * dyn3d_common: 1536 Bug correction concerning zoom (cf LMDZ5 rev 2218) 1537 - coefpoly.F becomes coefpoly_m.F90 (in misc) 1538 - fxhyp.F => fxhyp_m.F90 , fyhyp.F => fyhyp_m.F90 1539 - new routines for zoom: invert_zoom_x_m.F90 and principal_cshift_m.F90 1540 - inigeom.F adapted to new zoom definition routines 1541 - fluxstokenc.F : got rid of calls to initial0() 1542 1543 * dyn3d: 1544 - advtrac.F90 : got rid of calls to initial0() 1545 - conf_gcm.F90 : cosmetic changes and change in default dzoomx,dzoomy values 1546 - guide_mod.F90 : followed updates from Earth Model 1547 - gcm.F is now gcm.F90 1548 1549 * dyn3dpar: 1550 - advtrac_p.F90, covcont_p.F90, mod_hallo.F90 : cosmetic changes 1551 - conf_gcm.F90 : cosmetic and changed in default dzoomx,dzoomy values 1552 - parallel_lmdz.F90 : updates to keep up with Earth model 1553 1554 * misc: 1555 - arth.F90 becomes arth_m.F90 1556 - wxios.F90 updated wrt Earth model changes 1557 - nrtype.F90 and coefpoly_m.F90 added 1558 - ran1.F, sort.F, minmax.F, minmax2.F, juldate.F moved over from dyn3d_common 1559 -
trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90
r1422 r1441 74 74 75 75 IF(iadvtr.EQ.0) THEN 76 CALL initial0(ijp1llm,pbaruc)77 CALL initial0(ijmllm,pbarvc)76 pbaruc(:,:)=0 77 pbarvc(:,:)=0 78 78 ENDIF 79 79 -
trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F90
r1422 r1441 34 34 ! -metres du zoom avec celles lues sur le fichier start . 35 35 ! 36 LOGICAL etatinit37 INTEGER tapedef36 LOGICAL,INTENT(IN) :: etatinit 37 INTEGER,INTENT(IN) :: tapedef 38 38 39 39 ! Declarations : … … 44 44 include "iniprint.h" 45 45 46 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique47 ! #include "clesphys.h"48 !49 !50 46 ! local: 51 47 ! ------ … … 858 854 !Config Help = extension en longitude de la zone du zoom 859 855 !Config ( fraction de la zone totale) 860 dzoomx = 0. 0856 dzoomx = 0.2 861 857 CALL getin('dzoomx',dzoomx) 858 call assert(dzoomx < 1, "conf_gcm: dzoomx must be < 1") 862 859 863 860 !Config Key = dzoomy … … 866 863 !Config Help = extension en latitude de la zone du zoom 867 864 !Config ( fraction de la zone totale) 868 dzoomy = 0. 0865 dzoomy = 0.2 869 866 CALL getin('dzoomy',dzoomy) 867 call assert(dzoomy < 1, "conf_gcm: dzoomy must be < 1") 870 868 871 869 !Config Key = taux -
trunk/LMDZ.COMMON/libf/dyn3d/fluxstokenc.F
r1422 r1441 80 80 81 81 IF(iadvtr.EQ.0) THEN 82 CALL initial0(ijp1llm,phic)83 CALL initial0(ijp1llm,tetac)84 CALL initial0(ijp1llm,pbaruc)85 CALL initial0(ijmllm,pbarvc)82 phic(:,:)=0 83 tetac(:,:)=0 84 pbaruc(:,:)=0 85 pbarvc(:,:)=0 86 86 ENDIF 87 87 -
trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90
r1440 r1441 2 2 ! $Id: gcm.F 1446 2010-10-22 09:27:25Z emillour $ 3 3 ! 4 c 5 c 6 4 ! 5 ! 6 PROGRAM gcm 7 7 8 8 #ifdef CPP_IOIPSL 9 9 USE IOIPSL 10 10 #else 11 ! if not using IOIPSL, we still need to use (a local version of) getin12 11 ! if not using IOIPSL, we still need to use (a local version of) getin 12 USE ioipsl_getincom 13 13 #endif 14 14 15 15 16 16 #ifdef CPP_XIOS 17 18 19 #endif 20 21 22 23 USE control_mod, only: planet_type,nday,day_step,iperiod,iphysiq,24 & raz_date,anneeref,starttime,dayref,25 & ok_dyn_ins,ok_dyn_ave,iecri,periodav,26 &less1day,fractday,ndynstep,nsplit_phys27 28 USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref,29 .itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end17 ! ug Pour les sorties XIOS 18 USE wxios 19 #endif 20 21 USE filtreg_mod 22 USE infotrac 23 USE control_mod, only: planet_type,nday,day_step,iperiod,iphysiq, & 24 raz_date,anneeref,starttime,dayref, & 25 ok_dyn_ins,ok_dyn_ave,iecri,periodav, & 26 less1day,fractday,ndynstep,nsplit_phys 27 use cpdet_mod, only: ini_cpdet 28 USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref, & 29 itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end 30 30 31 31 #ifdef INCA 32 32 ! Only INCA needs these informations (from the Earth's physics) 33 USE indice_sol_mod 33 USE indice_sol_mod 34 USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb 34 35 #endif 35 36 … … 43 44 ! USE comgeomphy, ONLY: initcomgeomphy 44 45 #endif 45 #ifdef INCA 46 USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb 47 #endif 48 USE comconst_mod, ONLY: daysec,dtvr,dtphys,rad,g,r,cpp 49 USE logic_mod, ONLY: read_start,iflag_phys,ok_guide,ecripar 46 47 USE comconst_mod, ONLY: daysec,dtvr,dtphys,rad,g,r,cpp 48 USE logic_mod, ONLY: read_start,iflag_phys,ok_guide,ecripar 50 49 51 50 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 52 51 53 54 55 c...... Version du 10/01/98 ..........56 57 cavec coordonnees verticales hybrides58 cavec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )59 60 c=======================================================================61 c 62 cAuteur: P. Le Van /L. Fairhead/F.Hourdin63 c-------64 c 65 cObjet:66 c------67 c 68 cGCM LMD nouvelle grille69 c 70 c=======================================================================71 c 72 c... Dans inigeom , nouveaux calculs pour les elongations cu , cv73 cet possibilite d'appeler une fonction f(y) a derivee tangente74 chyperbolique a la place de la fonction a derivee sinusoidale.75 c... Possibilite de choisir le schema pour l'advection de76 cq , en modifiant iadv dans traceur.def (MAF,10/02) .77 c 78 cPour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)79 cPour Van-Leer iadv=1080 c 81 c-----------------------------------------------------------------------82 cDeclarations:83 c-------------84 85 #include "dimensions.h"86 #include "paramet.h"87 #include "comdissnew.h"88 #include "comgeom.h"52 IMPLICIT NONE 53 54 ! ...... Version du 10/01/98 .......... 55 56 ! avec coordonnees verticales hybrides 57 ! avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 ) 58 59 !======================================================================= 60 ! 61 ! Auteur: P. Le Van /L. Fairhead/F.Hourdin 62 ! ------- 63 ! 64 ! Objet: 65 ! ------ 66 ! 67 ! GCM LMD nouvelle grille 68 ! 69 !======================================================================= 70 ! 71 ! ... Dans inigeom , nouveaux calculs pour les elongations cu , cv 72 ! et possibilite d'appeler une fonction f(y) a derivee tangente 73 ! hyperbolique a la place de la fonction a derivee sinusoidale. 74 ! ... Possibilite de choisir le schema pour l'advection de 75 ! q , en modifiant iadv dans traceur.def (MAF,10/02) . 76 ! 77 ! Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99) 78 ! Pour Van-Leer iadv=10 79 ! 80 !----------------------------------------------------------------------- 81 ! Declarations: 82 ! ------------- 83 84 include "dimensions.h" 85 include "paramet.h" 86 include "comdissnew.h" 87 include "comgeom.h" 89 88 !!!!!!!!!!!#include "control.h" 90 89 !#include "com_io_dyn.h" 91 #include "iniprint.h"92 #include "tracstoke.h"90 include "iniprint.h" 91 include "tracstoke.h" 93 92 #ifdef INCA 94 93 ! Only INCA needs these informations (from the Earth's physics) … … 97 96 98 97 99 REAL zdtvr 100 101 c variables dynamiques 102 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 103 REAL teta(ip1jmp1,llm) ! temperature potentielle 104 REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes 105 REAL ps(ip1jmp1) ! pression au sol 106 REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches 107 REAL masse(ip1jmp1,llm) ! masse d'air 108 REAL phis(ip1jmp1) ! geopotentiel au sol 109 REAL phi(ip1jmp1,llm) ! geopotentiel 110 REAL w(ip1jmp1,llm) ! vitesse verticale 111 112 c variables dynamiques intermediaire pour le transport 113 114 c variables pour le fichier histoire 115 REAL dtav ! intervalle de temps elementaire 116 117 REAL time_0 118 119 LOGICAL lafin 120 INTEGER ij,iq,l,i,j 121 122 123 real time_step, t_wrt, t_ops 124 125 LOGICAL first 126 127 ! LOGICAL call_iniphys 128 ! data call_iniphys/.true./ 129 130 c+jld variables test conservation energie 131 c REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm) 132 C Tendance de la temp. potentiel d (theta)/ d t due a la 133 C tansformation d'energie cinetique en energie thermique 134 C cree par la dissipation 135 REAL dhecdt(ip1jmp1,llm) 136 c REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm) 137 c REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec 138 CHARACTER (len=15) :: ztit 139 c-jld 140 141 142 character (len=80) :: dynhist_file, dynhistave_file 143 character (len=20) :: modname 144 character (len=80) :: abort_message 145 ! locales pour gestion du temps 146 INTEGER :: an, mois, jour 147 REAL :: heure 148 149 150 c----------------------------------------------------------------------- 151 c variables pour l'initialisation de la physique : 152 c ------------------------------------------------ 153 ! INTEGER ngridmx 154 ! PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm ) 155 ! REAL zcufi(ngridmx),zcvfi(ngridmx) 156 ! REAL latfi(ngridmx),lonfi(ngridmx) 157 ! REAL airefi(ngridmx) 158 ! SAVE latfi, lonfi, airefi 159 160 c----------------------------------------------------------------------- 161 c Initialisations: 162 c ---------------- 163 164 abort_message = 'last timestep reached' 165 modname = 'gcm' 166 lafin = .FALSE. 167 dynhist_file = 'dyn_hist.nc' 168 dynhistave_file = 'dyn_hist_ave.nc' 169 170 171 172 c---------------------------------------------------------------------- 173 c lecture des fichiers gcm.def ou run.def 174 c --------------------------------------- 175 c 176 ! Ehouarn: dump possibility of using defrun 177 !#ifdef CPP_IOIPSL 178 CALL conf_gcm( 99, .TRUE. ) 179 if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", 180 s "iphysiq must be a multiple of iperiod", 1) 181 !#else 182 ! CALL defrun( 99, .TRUE. , clesphy0 ) 183 !#endif 98 REAL zdtvr 99 100 ! variables dynamiques 101 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 102 REAL teta(ip1jmp1,llm) ! temperature potentielle 103 REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes 104 REAL ps(ip1jmp1) ! pression au sol 105 REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches 106 REAL masse(ip1jmp1,llm) ! masse d'air 107 REAL phis(ip1jmp1) ! geopotentiel au sol 108 REAL phi(ip1jmp1,llm) ! geopotentiel 109 REAL w(ip1jmp1,llm) ! vitesse verticale 110 111 ! variables dynamiques intermediaire pour le transport 112 113 ! variables pour le fichier histoire 114 REAL dtav ! intervalle de temps elementaire 115 116 REAL time_0 117 118 LOGICAL lafin 119 INTEGER ij,iq,l,i,j 120 121 122 real time_step, t_wrt, t_ops 123 124 LOGICAL first 125 126 ! LOGICAL call_iniphys 127 ! data call_iniphys/.true./ 128 129 !+jld variables test conservation energie 130 ! REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm) 131 ! Tendance de la temp. potentiel d (theta)/ d t due a la 132 ! tansformation d'energie cinetique en energie thermique 133 ! cree par la dissipation 134 REAL dhecdt(ip1jmp1,llm) 135 ! REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm) 136 ! REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec 137 CHARACTER (len=15) :: ztit 138 !-jld 139 140 141 character (len=80) :: dynhist_file, dynhistave_file 142 character (len=20) :: modname 143 character (len=80) :: abort_message 144 ! locales pour gestion du temps 145 INTEGER :: an, mois, jour 146 REAL :: heure 147 logical use_filtre_fft 148 149 !----------------------------------------------------------------------- 150 ! Initialisations: 151 ! ---------------- 152 153 abort_message = 'last timestep reached' 154 modname = 'gcm' 155 lafin = .FALSE. 156 dynhist_file = 'dyn_hist.nc' 157 dynhistave_file = 'dyn_hist_ave.nc' 158 159 160 161 !---------------------------------------------------------------------- 162 ! lecture des fichiers gcm.def ou run.def 163 ! --------------------------------------- 164 ! 165 CALL conf_gcm( 99, .TRUE. ) 166 if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", & 167 "iphysiq must be a multiple of iperiod", 1) 168 169 use_filtre_fft=.FALSE. 170 CALL getin('use_filtre_fft',use_filtre_fft) 171 IF (use_filtre_fft) call abort_gcm('FFT filter is not available in the ' & 172 // 'sequential version of the dynamics.', 1) 184 173 185 174 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 188 177 189 178 #ifdef CPP_XIOS 190 179 CALL wxios_init("LMDZ") 191 180 #endif 192 181 … … 197 186 ! dynamique -> physique pour l'initialisation 198 187 #ifdef CPP_PHYS 199 188 CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) 200 189 ! call initcomgeomphy ! now done in iniphysiq 201 190 #endif 202 191 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 203 c 204 cInitialisations pour Cp(T) Venus205 206 c 207 c-----------------------------------------------------------------------208 cChoix du calendrier209 c-------------------210 211 ccalend = 'earth_365d'192 ! 193 ! Initialisations pour Cp(T) Venus 194 call ini_cpdet 195 ! 196 !----------------------------------------------------------------------- 197 ! Choix du calendrier 198 ! ------------------- 199 200 ! calend = 'earth_365d' 212 201 213 202 #ifdef CPP_IOIPSL 214 215 216 217 218 219 220 221 222 223 203 if (calend == 'earth_360d') then 204 call ioconf_calendar('360d') 205 write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an' 206 else if (calend == 'earth_365d') then 207 call ioconf_calendar('noleap') 208 write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an' 209 else if (calend == 'earth_366d') then 210 call ioconf_calendar('gregorian') 211 write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile' 212 else if (calend == 'titan') then 224 213 ! call ioconf_calendar('titan') 225 226 227 228 214 write(lunout,*)'CALENDRIER CHOISI: Titan' 215 abort_message = 'A FAIRE...' 216 call abort_gcm(modname,abort_message,1) 217 else if (calend == 'venus') then 229 218 ! call ioconf_calendar('venus') 230 231 232 233 234 235 236 237 #endif 238 c-----------------------------------------------------------------------239 240 219 write(lunout,*)'CALENDRIER CHOISI: Venus' 220 abort_message = 'A FAIRE...' 221 call abort_gcm(modname,abort_message,1) 222 else 223 abort_message = 'Mauvais choix de calendrier' 224 call abort_gcm(modname,abort_message,1) 225 endif 226 #endif 227 !----------------------------------------------------------------------- 228 229 IF (type_trac == 'inca') THEN 241 230 #ifdef INCA 242 call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,243 $nbsrf, is_oce,is_sic,is_ter,is_lic)244 245 #endif 246 247 c 248 c 249 c------------------------------------250 cInitialisation partie parallele251 c------------------------------------252 253 c 254 c 255 c-----------------------------------------------------------------------256 cInitialisation des traceurs257 c---------------------------258 cChoix du nombre de traceurs et du schema pour l'advection259 cdans fichier traceur.def, par default ou via INCA260 261 262 cAllocation de la tableau q : champs advectes263 264 265 c-----------------------------------------------------------------------266 cLecture de l'etat initial :267 c---------------------------268 269 clecture du fichier start.nc270 271 272 273 274 275 276 277 CALL dynetat0("start.nc",vcov,ucov,278 &teta,q,masse,ps,phis, time_0)231 call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday, & 232 nbsrf, is_oce,is_sic,is_ter,is_lic) 233 call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0) 234 #endif 235 END IF 236 ! 237 ! 238 !------------------------------------ 239 ! Initialisation partie parallele 240 !------------------------------------ 241 242 ! 243 ! 244 !----------------------------------------------------------------------- 245 ! Initialisation des traceurs 246 ! --------------------------- 247 ! Choix du nombre de traceurs et du schema pour l'advection 248 ! dans fichier traceur.def, par default ou via INCA 249 call infotrac_init 250 251 ! Allocation de la tableau q : champs advectes 252 allocate(q(ip1jmp1,llm,nqtot)) 253 254 !----------------------------------------------------------------------- 255 ! Lecture de l'etat initial : 256 ! --------------------------- 257 258 ! lecture du fichier start.nc 259 if (read_start) then 260 ! we still need to run iniacademic to initialize some 261 ! constants & fields, if we run the 'newtonian' or 'SW' cases: 262 if (iflag_phys.ne.1) then 263 CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) 264 endif 265 266 CALL dynetat0("start.nc",vcov,ucov, & 267 teta,q,masse,ps,phis, time_0) 279 268 280 281 282 283 284 285 286 269 ! Load relaxation fields (simple nudging). AS 09/2013 270 ! --------------------------------------------------- 271 if (planet_type.eq."generic") then 272 if (ok_guide) then 273 CALL relaxetat0("relax.nc") 274 endif 275 endif 287 276 288 cwrite(73,*) 'ucov',ucov289 cwrite(74,*) 'vcov',vcov290 cwrite(75,*) 'teta',teta291 cwrite(76,*) 'ps',ps292 cwrite(77,*) 'q',q293 294 295 296 277 ! write(73,*) 'ucov',ucov 278 ! write(74,*) 'vcov',vcov 279 ! write(75,*) 'teta',teta 280 ! write(76,*) 'ps',ps 281 ! write(77,*) 'q',q 282 283 endif ! of if (read_start) 284 285 IF (type_trac == 'inca') THEN 297 286 #ifdef INCA 298 call init_inca_dim(klon,llm,iim,jjm,299 $rlonu,rlatu,rlonv,rlatv)300 #endif 301 302 303 304 cle cas echeant, creation d un etat initial305 IF (prt_level > 9) WRITE(lunout,*)306 .'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'307 308 309 310 311 312 c-----------------------------------------------------------------------313 cLecture des parametres de controle pour la simulation :314 c-------------------------------------------------------315 con recalcule eventuellement le pas de temps316 317 318 abort_message =319 .'Il faut choisir un nb de pas par jour multiple de iperiod'320 321 322 323 324 abort_message =325 *'Il faut choisir un nb de pas par jour multiple de iphysiq'326 327 328 329 330 331 WRITE(lunout,*)332 .'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr333 334 335 C 336 Con remet le calendrier à zero si demande337 c 338 339 WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'340 &,' fichier restart ne correspond pas à celle lue dans le run.def'341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 write(lunout,*)356 .'GCM: On reinitialise a la date lue dans gcm.def'357 358 write(lunout,*)359 .'GCM: Attention les dates initiales lues dans le fichier'360 write(lunout,*)361 .' restart ne correspondent pas a celles lues dans '362 363 364 365 366 367 368 cif (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then369 cwrite(lunout,*)370 c. 'GCM: Attention les dates initiales lues dans le fichier'371 cwrite(lunout,*)372 c. ' restart ne correspondent pas a celles lues dans '373 cwrite(lunout,*)' gcm.def'374 cwrite(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref375 cwrite(lunout,*)' day_ref=',day_ref," dayref=",dayref376 cif (raz_date .ne. 1) then377 cwrite(lunout,*)378 c. 'GCM: On garde les dates du fichier restart'379 celse380 cannee_ref = anneeref381 cday_ref = dayref382 cday_ini = dayref383 citau_dyn = 0384 citau_phy = 0385 ctime_0 = 0.386 cwrite(lunout,*)387 c. 'GCM: On reinitialise a la date lue dans gcm.def'388 cendif389 cELSE390 craz_date = 0391 cendif287 call init_inca_dim(klon,llm,iim,jjm, & 288 rlonu,rlatu,rlonv,rlatv) 289 #endif 290 END IF 291 292 293 ! le cas echeant, creation d un etat initial 294 IF (prt_level > 9) WRITE(lunout,*) & 295 'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT' 296 if (.not.read_start) then 297 CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) 298 endif 299 300 301 !----------------------------------------------------------------------- 302 ! Lecture des parametres de controle pour la simulation : 303 ! ------------------------------------------------------- 304 ! on recalcule eventuellement le pas de temps 305 306 IF(MOD(day_step,iperiod).NE.0) THEN 307 abort_message = & 308 'Il faut choisir un nb de pas par jour multiple de iperiod' 309 call abort_gcm(modname,abort_message,1) 310 ENDIF 311 312 IF(MOD(day_step,iphysiq).NE.0) THEN 313 abort_message = & 314 'Il faut choisir un nb de pas par jour multiple de iphysiq' 315 call abort_gcm(modname,abort_message,1) 316 ENDIF 317 318 zdtvr = daysec/REAL(day_step) 319 IF(dtvr.NE.zdtvr) THEN 320 WRITE(lunout,*) & 321 'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr 322 ENDIF 323 324 ! 325 ! on remet le calendrier à zero si demande 326 ! 327 IF (start_time /= starttime) then 328 WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' & 329 ,' fichier restart ne correspond pas à celle lue dans le run.def' 330 IF (raz_date == 1) then 331 WRITE(lunout,*)'Je prends l''heure lue dans run.def' 332 start_time = starttime 333 ELSE 334 call abort_gcm("gcm", "'Je m''arrete'", 1) 335 ENDIF 336 ENDIF 337 IF (raz_date == 1) THEN 338 annee_ref = anneeref 339 day_ref = dayref 340 day_ini = dayref 341 itau_dyn = 0 342 itau_phy = 0 343 time_0 = 0. 344 write(lunout,*) & 345 'GCM: On reinitialise a la date lue dans gcm.def' 346 ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN 347 write(lunout,*) & 348 'GCM: Attention les dates initiales lues dans le fichier' 349 write(lunout,*) & 350 ' restart ne correspondent pas a celles lues dans ' 351 write(lunout,*)' gcm.def' 352 write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref 353 write(lunout,*)' day_ref=',day_ref," dayref=",dayref 354 write(lunout,*)' Pas de remise a zero' 355 ENDIF 356 357 ! if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then 358 ! write(lunout,*) 359 ! . 'GCM: Attention les dates initiales lues dans le fichier' 360 ! write(lunout,*) 361 ! . ' restart ne correspondent pas a celles lues dans ' 362 ! write(lunout,*)' gcm.def' 363 ! write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref 364 ! write(lunout,*)' day_ref=',day_ref," dayref=",dayref 365 ! if (raz_date .ne. 1) then 366 ! write(lunout,*) 367 ! . 'GCM: On garde les dates du fichier restart' 368 ! else 369 ! annee_ref = anneeref 370 ! day_ref = dayref 371 ! day_ini = dayref 372 ! itau_dyn = 0 373 ! itau_phy = 0 374 ! time_0 = 0. 375 ! write(lunout,*) 376 ! . 'GCM: On reinitialise a la date lue dans gcm.def' 377 ! endif 378 ! ELSE 379 ! raz_date = 0 380 ! endif 392 381 393 382 #ifdef CPP_IOIPSL 394 395 383 mois = 1 384 heure = 0. 396 385 ! Ce n'est defini pour l'instant que pour la Terre... 397 398 399 400 401 402 403 404 405 406 407 408 409 410 386 if (planet_type.eq.'earth') then 387 call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref) 388 jH_ref = jD_ref - int(jD_ref) 389 jD_ref = int(jD_ref) 390 391 call ioconf_startdate(INT(jD_ref), jH_ref) 392 393 write(lunout,*)'DEBUG' 394 write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref' 395 write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref 396 call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure) 397 write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure' 398 write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure 399 else 411 400 ! A voir pour Titan et Venus 412 413 414 415 416 401 jD_ref=0 402 jH_ref=0 403 write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref' 404 write(lunout,*)jD_ref,jH_ref 405 endif ! planet_type 417 406 #else 418 ! Ehouarn: we still need to define JD_ref and JH_ref 419 ! and since we don't know how many days there are in a year 420 ! we set JD_ref to 0 (this should be improved ...) 421 jD_ref=0 422 jH_ref=0 423 #endif 424 425 if (iflag_phys.eq.1) then 426 ! these initialisations have already been done (via iniacademic) 427 ! if running in SW or Newtonian mode 428 c----------------------------------------------------------------------- 429 c Initialisation des constantes dynamiques : 430 c ------------------------------------------ 431 dtvr = zdtvr 432 CALL iniconst 433 434 c----------------------------------------------------------------------- 435 c Initialisation de la geometrie : 436 c -------------------------------- 437 CALL inigeom 438 439 c----------------------------------------------------------------------- 440 c Initialisation du filtre : 441 c -------------------------- 442 CALL inifilr 443 endif ! of if (iflag_phys.eq.1) 444 c 445 c----------------------------------------------------------------------- 446 c Initialisation de la dissipation : 447 c ---------------------------------- 448 449 CALL inidissip( lstardis, nitergdiv, nitergrot, niterh , 450 * tetagdiv, tetagrot , tetatemp, vert_prof_dissip) 451 452 c----------------------------------------------------------------------- 453 c Initialisation de la physique : 454 c ------------------------------- 455 456 IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN 457 ! IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN 458 ! latfi(1)=rlatu(1) 459 ! lonfi(1)=0. 460 ! zcufi(1) = cu(1) 461 ! zcvfi(1) = cv(1) 462 ! DO j=2,jjm 463 ! DO i=1,iim 464 ! latfi((j-2)*iim+1+i)= rlatu(j) 465 ! lonfi((j-2)*iim+1+i)= rlonv(i) 466 ! zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i) 467 ! zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i) 468 ! ENDDO 469 ! ENDDO 470 ! latfi(ngridmx)= rlatu(jjp1) 471 ! lonfi(ngridmx)= 0. 472 ! zcufi(ngridmx) = cu(ip1jm+1) 473 ! zcvfi(ngridmx) = cv(ip1jm-iim) 474 475 ! build airefi(), mesh area on physics grid 476 ! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi) 477 ! Poles are single points on physics grid 478 ! airefi(1)=airefi(1)*iim 479 ! airefi(ngridmx)=airefi(ngridmx)*iim 480 481 ! Initialisation de la physique: pose probleme quand on tourne 482 ! SANS physique, car iniphysiq.F est dans le repertoire phy[]... 483 ! Il faut une cle CPP_PHYS 407 ! Ehouarn: we still need to define JD_ref and JH_ref 408 ! and since we don't know how many days there are in a year 409 ! we set JD_ref to 0 (this should be improved ...) 410 jD_ref=0 411 jH_ref=0 412 #endif 413 414 if (iflag_phys.eq.1) then 415 ! these initialisations have already been done (via iniacademic) 416 ! if running in SW or Newtonian mode 417 !----------------------------------------------------------------------- 418 ! Initialisation des constantes dynamiques : 419 ! ------------------------------------------ 420 dtvr = zdtvr 421 CALL iniconst 422 423 !----------------------------------------------------------------------- 424 ! Initialisation de la geometrie : 425 ! -------------------------------- 426 CALL inigeom 427 428 !----------------------------------------------------------------------- 429 ! Initialisation du filtre : 430 ! -------------------------- 431 CALL inifilr 432 endif ! of if (iflag_phys.eq.1) 433 ! 434 !----------------------------------------------------------------------- 435 ! Initialisation de la dissipation : 436 ! ---------------------------------- 437 438 CALL inidissip( lstardis, nitergdiv, nitergrot, niterh , & 439 tetagdiv, tetagrot , tetatemp, vert_prof_dissip) 440 441 !----------------------------------------------------------------------- 442 ! Initialisation de la physique : 443 ! ------------------------------- 444 445 IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN 446 ! Physics: 484 447 #ifdef CPP_PHYS 485 ! CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys, 486 ! & latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp, 487 ! & iflag_phys) 488 CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, 489 & rlatu,rlonv,aire,cu,cv,rad,g,r,cpp, 490 & iflag_phys) 491 #endif 492 ! call_iniphys=.false. 493 ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100)) 494 495 c numero de stockage pour les fichiers de redemarrage: 496 497 c----------------------------------------------------------------------- 498 c Initialisation des I/O : 499 c ------------------------ 500 501 if (nday>=0) then ! standard case 502 day_end=day_ini+nday 503 else ! special case when nday <0, run -nday dynamical steps 504 day_end=day_ini-nday/day_step 505 endif 506 if (less1day) then 507 day_end=day_ini+floor(time_0+fractday) 508 endif 509 if (ndynstep.gt.0) then 510 day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step)) 511 endif 448 CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, & 449 rlatu,rlonv,aire,cu,cv,rad,g,r,cpp, & 450 iflag_phys) 451 #endif 452 ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100)) 453 454 !----------------------------------------------------------------------- 455 ! Initialisation des I/O : 456 ! ------------------------ 457 458 if (nday>=0) then ! standard case 459 day_end=day_ini+nday 460 else ! special case when nday <0, run -nday dynamical steps 461 day_end=day_ini-nday/day_step 462 endif 463 if (less1day) then 464 day_end=day_ini+floor(time_0+fractday) 465 endif 466 if (ndynstep.gt.0) then 467 day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step)) 468 endif 512 469 513 WRITE(lunout,'(a,i7,a,i7)')514 &"run from day ",day_ini," to day",day_end470 WRITE(lunout,'(a,i7,a,i7)') & 471 "run from day ",day_ini," to day",day_end 515 472 516 473 #ifdef CPP_IOIPSL 517 ! Ce n'est defini pour l'instant que pour la Terre... 518 if (planet_type.eq.'earth') then 519 call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure) 520 write (lunout,301)jour, mois, an 521 call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure) 522 write (lunout,302)jour, mois, an 523 else 524 ! A voir pour Titan et Venus 525 write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...' 526 endif ! planet_type 527 528 301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4) 529 302 FORMAT('1'/,15x,' au ', i2,'/',i2,'/',i4) 530 #endif 531 532 if (planet_type=="mars") then 533 ! For Mars we transmit day_ini 534 CALL dynredem0("restart.nc", day_ini, phis) 535 else 536 CALL dynredem0("restart.nc", day_end, phis) 537 endif 538 ecripar = .TRUE. 474 ! Ce n'est defini pour l'instant que pour la Terre... 475 if (planet_type.eq.'earth') then 476 call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure) 477 write (lunout,301)jour, mois, an 478 call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure) 479 write (lunout,302)jour, mois, an 480 else 481 ! A voir pour Titan et Venus 482 write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...' 483 endif ! planet_type 484 301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4) 485 302 FORMAT('1'/,15x,' au ', i2,'/',i2,'/',i4) 486 #endif 487 488 if (planet_type=="mars") then 489 ! For Mars we transmit day_ini 490 CALL dynredem0("restart.nc", day_ini, phis) 491 else 492 CALL dynredem0("restart.nc", day_end, phis) 493 endif 494 ecripar = .TRUE. 539 495 540 496 #ifdef CPP_IOIPSL 541 542 543 544 545 546 547 CALL inithist(day_ref,annee_ref,time_step,548 &t_ops,t_wrt)549 550 551 552 553 554 555 CALL initdynav(day_ref,annee_ref,time_step,556 &t_ops,t_wrt)557 558 497 time_step = zdtvr 498 if (ok_dyn_ins) then 499 ! initialize output file for instantaneous outputs 500 ! t_ops = iecri * daysec ! do operations every t_ops 501 t_ops =((1.0*iecri)/day_step) * daysec 502 t_wrt = daysec ! iecri * daysec ! write output every t_wrt 503 CALL inithist(day_ref,annee_ref,time_step, & 504 t_ops,t_wrt) 505 endif 506 507 IF (ok_dyn_ave) THEN 508 ! initialize output file for averaged outputs 509 t_ops = iperiod * time_step ! do operations every t_ops 510 t_wrt = periodav * daysec ! write output every t_wrt 511 CALL initdynav(day_ref,annee_ref,time_step, & 512 t_ops,t_wrt) 513 END IF 514 dtav = iperiod*dtvr/daysec 559 515 #endif 560 516 ! #endif of #ifdef CPP_IOIPSL 561 517 562 c Choix des frequences de stokage pour le offline 563 c istdyn=day_step/4 ! stockage toutes les 6h=1jour/4 564 c istdyn=day_step/12 ! stockage toutes les 2h=1jour/12 565 istdyn=day_step/4 ! stockage toutes les 6h=1jour/12 566 istphy=istdyn/iphysiq 567 568 569 c 570 c----------------------------------------------------------------------- 571 c Integration temporelle du modele : 572 c ---------------------------------- 573 574 c write(78,*) 'ucov',ucov 575 c write(78,*) 'vcov',vcov 576 c write(78,*) 'teta',teta 577 c write(78,*) 'ps',ps 578 c write(78,*) 'q',q 579 580 581 CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q, 582 . time_0) 583 584 END 585 518 ! Choix des frequences de stokage pour le offline 519 ! istdyn=day_step/4 ! stockage toutes les 6h=1jour/4 520 ! istdyn=day_step/12 ! stockage toutes les 2h=1jour/12 521 istdyn=day_step/4 ! stockage toutes les 6h=1jour/12 522 istphy=istdyn/iphysiq 523 524 525 ! 526 !----------------------------------------------------------------------- 527 ! Integration temporelle du modele : 528 ! ---------------------------------- 529 530 ! write(78,*) 'ucov',ucov 531 ! write(78,*) 'vcov',vcov 532 ! write(78,*) 'teta',teta 533 ! write(78,*) 'ps',ps 534 ! write(78,*) 'q',q 535 536 537 CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0) 538 539 END PROGRAM gcm 540 -
trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90
r1422 r1441 154 154 rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl) 155 155 if (rcod.NE.NF_NOERR) THEN 156 print *,'Guide: probleme -> pas de fichier apbp.nc'157 CALL abort_gcm(modname,abort_message,1)156 CALL abort_gcm(modname, & 157 'Guide: probleme -> pas de fichier apbp.nc',1) 158 158 endif 159 159 endif … … 163 163 rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl) 164 164 if (rcod.NE.NF_NOERR) THEN 165 print *,'Guide: probleme -> pas de fichier u.nc'166 CALL abort_gcm(modname,abort_message,1)165 CALL abort_gcm(modname, & 166 'Guide: probleme -> pas de fichier u.nc',1) 167 167 endif 168 168 endif … … 171 171 rcod=nf90_open('v.nc',nf90_nowrite,ncidpl) 172 172 if (rcod.NE.NF_NOERR) THEN 173 print *,'Guide: probleme -> pas de fichier v.nc'174 CALL abort_gcm(modname,abort_message,1)173 CALL abort_gcm(modname, & 174 'Guide: probleme -> pas de fichier v.nc',1) 175 175 endif 176 176 endif … … 179 179 rcod=nf90_open('T.nc',nf90_nowrite,ncidpl) 180 180 if (rcod.NE.NF_NOERR) THEN 181 print *,'Guide: probleme -> pas de fichier T.nc'182 CALL abort_gcm(modname,abort_message,1)181 CALL abort_gcm(modname, & 182 'Guide: probleme -> pas de fichier T.nc',1) 183 183 endif 184 184 endif … … 187 187 rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl) 188 188 if (rcod.NE.NF_NOERR) THEN 189 print *,'Guide: probleme -> pas de fichier hur.nc'190 CALL abort_gcm(modname,abort_message,1)189 CALL abort_gcm(modname, & 190 'Guide: probleme -> pas de fichier hur.nc',1) 191 191 endif 192 192 endif … … 196 196 IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid) 197 197 IF (error.NE.NF_NOERR) THEN 198 print *,'Guide: probleme lecture niveaux pression' 199 CALL abort_gcm(modname,abort_message,1) 198 CALL abort_gcm(modname,'Guide: probleme lecture niveaux pression',1) 200 199 ENDIF 201 200 error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc) -
trunk/LMDZ.COMMON/libf/dyn3d_common/fxhyp_m.F90
r1440 r1441 1 ! 2 ! $Id: fxhyp.F 1674 2012-10-29 16:27:03Z emillour $ 3 ! 4 c 5 c 6 SUBROUTINE fxhyp ( xzoomdeg,grossism,dzooma,tau , 7 , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025, 8 , champmin,champmax ) 9 10 c Auteur : P. Le Van 11 12 IMPLICIT NONE 13 14 c Calcule les longitudes et derivees dans la grille du GCM pour une 15 c fonction f(x) a tangente hyperbolique . 16 c 17 c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.) 18 c dzoom etant la distance totale de la zone du zoom 19 c tau la raideur de la transition de l'interieur a l'exterieur du zoom 20 c 21 c On doit avoir grossism x dzoom < pi ( radians ) , en longitude. 22 c ******************************************************************** 23 24 25 INTEGER nmax, nmax2 26 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 27 c 28 LOGICAL scal180 29 PARAMETER ( scal180 = .TRUE. ) 30 31 c scal180 = .TRUE. si on veut avoir le premier point scalaire pour 32 c une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres. 33 c sinon scal180 = .FALSE. 34 35 #include "dimensions.h" 36 #include "paramet.h" 37 38 c ...... arguments d'entree ....... 39 c 40 REAL xzoomdeg,dzooma,tau,grossism 41 42 c ...... arguments de sortie ...... 43 44 REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1), 45 , rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1) 46 47 c .... variables locales .... 48 c 49 REAL dzoom 50 REAL(KIND=8) xlon(iip1),xprimm(iip1),xuv 51 REAL(KIND=8) xtild(0:nmax2) 52 REAL(KIND=8) fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2) 53 REAL(KIND=8) Xf(0:nmax2),xxpr(0:nmax2) 54 REAL(KIND=8) xvrai(iip1),xxprim(iip1) 55 REAL(KIND=8) pi,depi,epsilon,xzoom,fa,fb 56 REAL(KIND=8) Xf1, Xfi , a0,a1,a2,a3,xi2 57 INTEGER i,it,ik,iter,ii,idif,ii1,ii2 58 REAL(KIND=8) xi,xo1,xmoy,xlon2,fxm,Xprimin 59 REAL(KIND=8) champmin,champmax,decalx 60 INTEGER is2 61 SAVE is2 62 63 REAL(KIND=8) heavyside 64 65 pi = 2. * ASIN(1.) 66 depi = 2. * pi 67 epsilon = 1.e-3 68 xzoom = xzoomdeg * pi/180. 69 c 70 if (iim==1) then 71 72 rlonm025(1)=-pi/2. 73 rlonv(1)=0. 74 rlonu(1)=pi 75 rlonp025(1)=pi/2. 76 rlonm025(2)=rlonm025(1)+depi 77 rlonv(2)=rlonv(1)+depi 78 rlonu(2)=rlonu(1)+depi 79 rlonp025(2)=rlonp025(1)+depi 80 81 xprimm025(:)=1. 82 xprimv(:)=1. 83 xprimu(:)=1. 84 xprimp025(:)=1. 85 champmin=depi 86 champmax=depi 87 return 88 89 endif 90 91 decalx = .75 92 IF( grossism.EQ.1..AND.scal180 ) THEN 93 decalx = 1. 94 ENDIF 95 96 WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx 97 c 98 IF( dzooma.LT.1.) THEN 99 dzoom = dzooma * depi 100 ELSEIF( dzooma.LT. 25. ) THEN 101 WRITE(6,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug 102 ,menter et relancer ! ' 103 STOP 1 104 ELSE 105 dzoom = dzooma * pi/180. 106 ENDIF 107 108 WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)' 109 WRITE(6,24) xzoom,grossism,tau,dzoom 110 111 DO i = 0, nmax2 112 xtild(i) = - pi + REAL(i) * depi /nmax2 113 ENDDO 114 115 DO i = nmax, nmax2 116 117 fa = tau* ( dzoom/2. - xtild(i) ) 118 fb = xtild(i) * ( pi - xtild(i) ) 119 120 IF( 200.* fb .LT. - fa ) THEN 121 fhyp ( i) = - 1. 122 ELSEIF( 200. * fb .LT. fa ) THEN 123 fhyp ( i) = 1. 124 ELSE 125 IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN 126 IF( 200.*fb + fa.LT.1.e-10 ) THEN 127 fhyp ( i ) = - 1. 128 ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN 129 fhyp ( i ) = 1. 130 ENDIF 131 ELSE 132 fhyp ( i ) = TANH ( fa/fb ) 133 ENDIF 134 ENDIF 135 IF ( xtild(i).EQ. 0. ) fhyp(i) = 1. 136 IF ( xtild(i).EQ. pi ) fhyp(i) = -1. 137 138 ENDDO 139 140 cc .... Calcul de beta .... 141 142 ffdx = 0. 143 144 DO i = nmax +1,nmax2 145 146 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 147 fa = tau* ( dzoom/2. - xmoy ) 148 fb = xmoy * ( pi - xmoy ) 149 150 IF( 200.* fb .LT. - fa ) THEN 151 fxm = - 1. 152 ELSEIF( 200. * fb .LT. fa ) THEN 153 fxm = 1. 154 ELSE 155 IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN 156 IF( 200.*fb + fa.LT.1.e-10 ) THEN 157 fxm = - 1. 158 ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN 159 fxm = 1. 160 ENDIF 161 ELSE 162 fxm = TANH ( fa/fb ) 163 ENDIF 164 ENDIF 165 166 IF ( xmoy.EQ. 0. ) fxm = 1. 167 IF ( xmoy.EQ. pi ) fxm = -1. 168 169 ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) ) 170 171 ENDDO 172 173 beta = ( grossism * ffdx - pi ) / ( ffdx - pi ) 174 175 IF( 2.*beta - grossism.LE. 0.) THEN 176 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 177 ,tine fxhyp est mauvaise ! ' 178 WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ', 179 , ' et relancer ! *** ' 180 CALL ABORT_GCM("FXHYP", "", 1) 181 ENDIF 182 c 183 c ..... calcul de Xprimt ..... 184 c 185 186 DO i = nmax, nmax2 187 Xprimt(i) = beta + ( grossism - beta ) * fhyp(i) 188 ENDDO 189 c 190 DO i = nmax+1, nmax2 191 Xprimt( nmax2 - i ) = Xprimt( i ) 192 ENDDO 193 c 194 195 c ..... Calcul de Xf ........ 196 197 Xf(0) = - pi 198 199 DO i = nmax +1, nmax2 200 201 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 202 fa = tau* ( dzoom/2. - xmoy ) 203 fb = xmoy * ( pi - xmoy ) 204 205 IF( 200.* fb .LT. - fa ) THEN 206 fxm = - 1. 207 ELSEIF( 200. * fb .LT. fa ) THEN 208 fxm = 1. 209 ELSE 210 fxm = TANH ( fa/fb ) 211 ENDIF 212 213 IF ( xmoy.EQ. 0. ) fxm = 1. 214 IF ( xmoy.EQ. pi ) fxm = -1. 215 xxpr(i) = beta + ( grossism - beta ) * fxm 216 217 ENDDO 218 219 DO i = nmax+1, nmax2 220 xxpr(nmax2-i+1) = xxpr(i) 221 ENDDO 222 223 DO i=1,nmax2 224 Xf(i) = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) ) 225 ENDDO 226 227 228 c ***************************************************************** 229 c 230 231 c ..... xuv = 0. si calcul aux pts scalaires ........ 232 c ..... xuv = 0.5 si calcul aux pts U ........ 233 c 234 WRITE(6,18) 235 c 236 DO 5000 ik = 1, 4 237 238 IF( ik.EQ.1 ) THEN 239 xuv = -0.25 240 ELSE IF ( ik.EQ.2 ) THEN 241 xuv = 0. 242 ELSE IF ( ik.EQ.3 ) THEN 243 xuv = 0.50 244 ELSE IF ( ik.EQ.4 ) THEN 245 xuv = 0.25 246 ENDIF 247 248 xo1 = 0. 249 250 ii1=1 251 ii2=iim 252 IF(ik.EQ.1.and.grossism.EQ.1.) THEN 253 ii1 = 2 254 ii2 = iim+1 255 ENDIF 256 DO 1500 i = ii1, ii2 257 258 xlon2 = - pi + (REAL(i) + xuv - decalx) * depi / REAL(iim) 259 260 Xfi = xlon2 261 c 262 DO 250 it = nmax2,0,-1 263 IF( Xfi.GE.Xf(it)) GO TO 350 264 250 CONTINUE 265 266 it = 0 267 268 350 CONTINUE 269 270 c ...... Calcul de Xf(xi) ...... 271 c 272 xi = xtild(it) 273 274 IF(it.EQ.nmax2) THEN 275 it = nmax2 -1 276 Xf(it+1) = pi 277 ENDIF 278 c ..................................................................... 279 c 280 c Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un 281 c polynome de degre 3 qui passe par les points (Xf(it),xtild(it) ) 282 c et (Xf(it+1),xtild(it+1) ) 283 284 CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1), 285 , xtild(it),xtild(it+1), a0, a1, a2, a3 ) 286 287 Xf1 = Xf(it) 288 Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi 289 290 DO 500 iter = 1,300 291 xi = xi - ( Xf1 - Xfi )/ Xprimin 292 293 IF( ABS(xi-xo1).LE.epsilon) GO TO 550 294 xo1 = xi 295 xi2 = xi * xi 296 Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi 297 Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2 298 500 CONTINUE 299 WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter 300 STOP 6 301 550 CONTINUE 302 303 xxprim(i) = depi/ ( REAL(iim) * Xprimin ) 304 xvrai(i) = xi + xzoom 305 306 1500 CONTINUE 307 308 309 310 IF(ik.EQ.1.and.grossism.EQ.1.) THEN 311 xvrai(1) = xvrai(iip1)-depi 312 xxprim(1) = xxprim(iip1) 313 ENDIF 314 315 DO i = 1 , iim 316 xlon(i) = xvrai(i) 317 xprimm(i) = xxprim(i) 318 ENDDO 319 DO i = 1, iim -1 320 IF( xvrai(i+1). LT. xvrai(i) ) THEN 321 WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i, 322 , ')' 323 STOP 7 324 ENDIF 325 ENDDO 326 c 327 c ... Reorganisation des longitudes pour les avoir entre - pi et pi .. 328 c ........................................................................ 329 330 champmin = 1.e12 331 champmax = -1.e12 332 DO i = 1, iim 333 champmin = MIN( champmin,xvrai(i) ) 334 champmax = MAX( champmax,xvrai(i) ) 335 ENDDO 336 337 IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 ) THEN 338 GO TO 1600 339 ELSE 340 WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi', 341 , ' et pi ' 342 c 343 IF( xzoom.LE.0.) THEN 344 IF( ik.EQ. 1 ) THEN 345 DO i = 1, iim 346 IF( xvrai(i).GE. - pi ) GO TO 80 347 ENDDO 348 WRITE(6,*) ' PBS. 1 ! Xvrai plus petit que - pi ! ' 349 STOP 8 350 80 CONTINUE 351 is2 = i 352 ENDIF 353 354 IF( is2.NE. 1 ) THEN 355 DO ii = is2 , iim 356 xlon (ii-is2+1) = xvrai(ii) 357 xprimm(ii-is2+1) = xxprim(ii) 358 ENDDO 359 DO ii = 1 , is2 -1 360 xlon (ii+iim-is2+1) = xvrai(ii) + depi 361 xprimm(ii+iim-is2+1) = xxprim(ii) 362 ENDDO 363 ENDIF 364 ELSE 365 IF( ik.EQ.1 ) THEN 366 DO i = iim,1,-1 367 IF( xvrai(i).LE. pi ) GO TO 90 368 ENDDO 369 WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! ' 370 STOP 9 371 90 CONTINUE 372 is2 = i 373 ENDIF 374 idif = iim -is2 375 DO ii = 1, is2 376 xlon (ii+idif) = xvrai(ii) 377 xprimm(ii+idif) = xxprim(ii) 378 ENDDO 379 DO ii = 1, idif 380 xlon (ii) = xvrai (ii+is2) - depi 381 xprimm(ii) = xxprim(ii+is2) 382 ENDDO 383 ENDIF 384 ENDIF 385 c 386 c ......... Fin de la reorganisation ............................ 387 388 1600 CONTINUE 389 390 391 xlon ( iip1) = xlon(1) + depi 392 xprimm( iip1 ) = xprimm (1 ) 393 394 DO i = 1, iim+1 395 xvrai(i) = xlon(i)*180./pi 396 ENDDO 397 398 IF( ik.EQ.1 ) THEN 399 c WRITE(6,*) ' XLON aux pts. V-0.25 apres ( en deg. ) ' 400 c WRITE(6,18) 401 c WRITE(6,68) xvrai 402 c WRITE(6,*) ' XPRIM k ',ik 403 c WRITE(6,566) xprimm 404 405 DO i = 1,iim +1 406 rlonm025(i) = xlon( i ) 407 xprimm025(i) = xprimm(i) 408 ENDDO 409 ELSE IF( ik.EQ.2 ) THEN 410 c WRITE(6,18) 411 c WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) ' 412 c WRITE(6,68) xvrai 413 c WRITE(6,*) ' XPRIM k ',ik 414 c WRITE(6,566) xprimm 415 416 DO i = 1,iim + 1 417 rlonv(i) = xlon( i ) 418 xprimv(i) = xprimm(i) 419 ENDDO 420 421 ELSE IF( ik.EQ.3) THEN 422 c WRITE(6,18) 423 c WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) ' 424 c WRITE(6,68) xvrai 425 c WRITE(6,*) ' XPRIM ik ',ik 426 c WRITE(6,566) xprimm 427 428 DO i = 1,iim + 1 429 rlonu(i) = xlon( i ) 430 xprimu(i) = xprimm(i) 431 ENDDO 432 433 ELSE IF( ik.EQ.4 ) THEN 434 c WRITE(6,18) 435 c WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) ' 436 c WRITE(6,68) xvrai 437 c WRITE(6,*) ' XPRIM ik ',ik 438 c WRITE(6,566) xprimm 439 440 DO i = 1,iim + 1 441 rlonp025(i) = xlon( i ) 442 xprimp025(i) = xprimm(i) 443 ENDDO 444 445 ENDIF 446 447 5000 CONTINUE 448 c 449 WRITE(6,18) 450 c 451 c ........... fin de la boucle do 5000 ............ 452 453 DO i = 1, iim 454 xlon(i) = rlonv(i+1) - rlonv(i) 455 ENDDO 456 champmin = 1.e12 457 champmax = -1.e12 458 DO i = 1, iim 459 champmin = MIN( champmin, xlon(i) ) 460 champmax = MAX( champmax, xlon(i) ) 461 ENDDO 462 champmin = champmin * 180./pi 463 champmax = champmax * 180./pi 464 465 18 FORMAT(/) 466 24 FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3) 467 68 FORMAT(1x,7f9.2) 468 566 FORMAT(1x,7f9.4) 469 470 RETURN 471 END 1 module fxhyp_m 2 3 IMPLICIT NONE 4 5 contains 6 7 SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) 8 9 ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32 10 ! Author: P. Le Van, from formulas by R. Sadourny 11 12 ! Calcule les longitudes et dérivées dans la grille du GCM pour 13 ! une fonction f(x) à dérivée tangente hyperbolique. 14 15 ! Il vaut mieux avoir : grossismx \times dzoom < pi 16 17 ! Le premier point scalaire pour une grille regulière (grossismx = 18 ! 1., taux=0., clon=0.) est à - 180 degrés. 19 20 use arth_m, only: arth 21 use invert_zoom_x_m, only: invert_zoom_x, nmax 22 use nrtype, only: pi, pi_d, twopi, twopi_d, k8 23 use principal_cshift_m, only: principal_cshift 24 25 include "dimensions.h" 26 ! for iim 27 28 include "serre.h" 29 ! for clon, grossismx, dzoomx, taux 30 31 REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1) 32 real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1) 33 34 ! Local: 35 real rlonm025(iim + 1), rlonp025(iim + 1) 36 REAL dzoom, step 37 real d_rlonv(iim) 38 REAL(K8) xtild(0:2 * nmax) 39 REAL(K8) fhyp(nmax:2 * nmax), ffdx, beta, Xprimt(0:2 * nmax) 40 REAL(K8) Xf(0:2 * nmax), xxpr(2 * nmax) 41 REAL(K8) fa, fb 42 INTEGER i, is2 43 REAL(K8) xmoy, fxm 44 45 !---------------------------------------------------------------------- 46 47 print *, "Call sequence information: fxhyp" 48 49 test_iim: if (iim==1) then 50 rlonv(1)=0. 51 rlonu(1)=pi 52 rlonv(2)=rlonv(1)+twopi 53 rlonu(2)=rlonu(1)+twopi 54 55 xprimm025(:)=1. 56 xprimv(:)=1. 57 xprimu(:)=1. 58 xprimp025(:)=1. 59 else test_iim 60 test_grossismx: if (grossismx == 1.) then 61 step = twopi / iim 62 63 xprimm025(:iim) = step 64 xprimp025(:iim) = step 65 xprimv(:iim) = step 66 xprimu(:iim) = step 67 68 rlonv(:iim) = arth(- pi + clon / 180. * pi, step, iim) 69 rlonm025(:iim) = rlonv(:iim) - 0.25 * step 70 rlonp025(:iim) = rlonv(:iim) + 0.25 * step 71 rlonu(:iim) = rlonv(:iim) + 0.5 * step 72 else test_grossismx 73 dzoom = dzoomx * twopi_d 74 xtild = arth(- pi_d, pi_d / nmax, 2 * nmax + 1) 75 76 ! Compute fhyp: 77 DO i = nmax, 2 * nmax 78 fa = taux * (dzoom / 2. - xtild(i)) 79 fb = xtild(i) * (pi_d - xtild(i)) 80 81 IF (200. * fb < - fa) THEN 82 fhyp(i) = - 1. 83 ELSE IF (200. * fb < fa) THEN 84 fhyp(i) = 1. 85 ELSE 86 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN 87 IF (200. * fb + fa < 1e-10) THEN 88 fhyp(i) = - 1. 89 ELSE IF (200. * fb - fa < 1e-10) THEN 90 fhyp(i) = 1. 91 END IF 92 ELSE 93 fhyp(i) = TANH(fa / fb) 94 END IF 95 END IF 96 97 IF (xtild(i) == 0.) fhyp(i) = 1. 98 IF (xtild(i) == pi_d) fhyp(i) = -1. 99 END DO 100 101 ! Calcul de beta 102 103 ffdx = 0. 104 105 DO i = nmax + 1, 2 * nmax 106 xmoy = 0.5 * (xtild(i-1) + xtild(i)) 107 fa = taux * (dzoom / 2. - xmoy) 108 fb = xmoy * (pi_d - xmoy) 109 110 IF (200. * fb < - fa) THEN 111 fxm = - 1. 112 ELSE IF (200. * fb < fa) THEN 113 fxm = 1. 114 ELSE 115 IF (ABS(fa) < 1e-13 .AND. ABS(fb) < 1e-13) THEN 116 IF (200. * fb + fa < 1e-10) THEN 117 fxm = - 1. 118 ELSE IF (200. * fb - fa < 1e-10) THEN 119 fxm = 1. 120 END IF 121 ELSE 122 fxm = TANH(fa / fb) 123 END IF 124 END IF 125 126 IF (xmoy == 0.) fxm = 1. 127 IF (xmoy == pi_d) fxm = -1. 128 129 ffdx = ffdx + fxm * (xtild(i) - xtild(i-1)) 130 END DO 131 132 print *, "ffdx = ", ffdx 133 beta = (grossismx * ffdx - pi_d) / (ffdx - pi_d) 134 print *, "beta = ", beta 135 136 IF (2. * beta - grossismx <= 0.) THEN 137 print *, 'Bad choice of grossismx, taux, dzoomx.' 138 print *, 'Decrease dzoomx or grossismx.' 139 STOP 1 140 END IF 141 142 ! calcul de Xprimt 143 Xprimt(nmax:2 * nmax) = beta + (grossismx - beta) * fhyp 144 xprimt(:nmax - 1) = xprimt(2 * nmax:nmax + 1:- 1) 145 146 ! Calcul de Xf 147 148 DO i = nmax + 1, 2 * nmax 149 xmoy = 0.5 * (xtild(i-1) + xtild(i)) 150 fa = taux * (dzoom / 2. - xmoy) 151 fb = xmoy * (pi_d - xmoy) 152 153 IF (200. * fb < - fa) THEN 154 fxm = - 1. 155 ELSE IF (200. * fb < fa) THEN 156 fxm = 1. 157 ELSE 158 fxm = TANH(fa / fb) 159 END IF 160 161 IF (xmoy == 0.) fxm = 1. 162 IF (xmoy == pi_d) fxm = -1. 163 xxpr(i) = beta + (grossismx - beta) * fxm 164 END DO 165 166 xxpr(:nmax) = xxpr(2 * nmax:nmax + 1:- 1) 167 168 Xf(0) = - pi_d 169 170 DO i=1, 2 * nmax - 1 171 Xf(i) = Xf(i-1) + xxpr(i) * (xtild(i) - xtild(i-1)) 172 END DO 173 174 Xf(2 * nmax) = pi_d 175 176 call invert_zoom_x(xf, xtild, Xprimt, rlonm025(:iim), & 177 xprimm025(:iim), xuv = - 0.25_k8) 178 call invert_zoom_x(xf, xtild, Xprimt, rlonv(:iim), xprimv(:iim), & 179 xuv = 0._k8) 180 call invert_zoom_x(xf, xtild, Xprimt, rlonu(:iim), xprimu(:iim), & 181 xuv = 0.5_k8) 182 call invert_zoom_x(xf, xtild, Xprimt, rlonp025(:iim), & 183 xprimp025(:iim), xuv = 0.25_k8) 184 end if test_grossismx 185 186 is2 = 0 187 188 IF (MINval(rlonm025(:iim)) < - pi - 0.1 & 189 .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN 190 IF (clon <= 0.) THEN 191 is2 = 1 192 193 do while (rlonm025(is2) < - pi .and. is2 < iim) 194 is2 = is2 + 1 195 end do 196 197 if (rlonm025(is2) < - pi) then 198 print *, 'Rlonm025 plus petit que - pi !' 199 STOP 1 200 end if 201 ELSE 202 is2 = iim 203 204 do while (rlonm025(is2) > pi .and. is2 > 1) 205 is2 = is2 - 1 206 end do 207 208 if (rlonm025(is2) > pi) then 209 print *, 'Rlonm025 plus grand que pi !' 210 STOP 1 211 end if 212 END IF 213 END IF 214 215 call principal_cshift(is2, rlonm025, xprimm025) 216 call principal_cshift(is2, rlonv, xprimv) 217 call principal_cshift(is2, rlonu, xprimu) 218 call principal_cshift(is2, rlonp025, xprimp025) 219 220 forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i) 221 print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, & 222 "degrees" 223 print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, & 224 "degrees" 225 226 ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu: 227 DO i = 1, iim + 1 228 IF (rlonp025(i) < rlonv(i)) THEN 229 print *, 'rlonp025(', i, ') = ', rlonp025(i) 230 print *, "< rlonv(", i, ") = ", rlonv(i) 231 STOP 1 232 END IF 233 234 IF (rlonv(i) < rlonm025(i)) THEN 235 print *, 'rlonv(', i, ') = ', rlonv(i) 236 print *, "< rlonm025(", i, ") = ", rlonm025(i) 237 STOP 1 238 END IF 239 240 IF (rlonp025(i) > rlonu(i)) THEN 241 print *, 'rlonp025(', i, ') = ', rlonp025(i) 242 print *, "> rlonu(", i, ") = ", rlonu(i) 243 STOP 1 244 END IF 245 END DO 246 end if test_iim 247 248 END SUBROUTINE fxhyp 249 250 end module fxhyp_m -
trunk/LMDZ.COMMON/libf/dyn3d_common/fyhyp_m.F90
r1440 r1441 1 ! 2 ! $Id: fyhyp.F 1403 2010-07-01 09:02:53Z fairhead $ 3 ! 4 c 5 c 6 SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau , 7 , rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 , 8 , champmin,champmax ) 9 10 cc ... Version du 01/04/2001 .... 11 12 IMPLICIT NONE 13 c 14 c ... Auteur : P. Le Van ... 15 c 16 c ....... d'apres formulations de R. Sadourny ....... 17 c 18 c Calcule les latitudes et derivees dans la grille du GCM pour une 19 c fonction f(y) a tangente hyperbolique . 20 c 21 c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc) 22 c dzoom etant la distance totale de la zone du zoom ( en radians ) 23 c tau la raideur de la transition de l'interieur a l'exterieur du zoom 24 c 25 c 26 c N.B : Il vaut mieux avoir : grossism * dzoom < pi/2 (radians) ,en lati. 27 c ******************************************************************** 28 c 29 c 30 #include "dimensions.h" 31 #include "paramet.h" 32 33 INTEGER nmax , nmax2 34 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 35 c 36 c 37 c ....... arguments d'entree ....... 38 c 39 REAL yzoomdeg, grossism,dzooma,tau 40 c ( rentres par run.def ) 41 42 c ....... arguments de sortie ....... 43 c 44 REAL rrlatu(jjp1), yyprimu(jjp1),rrlatv(jjm), yyprimv(jjm), 45 , rlatu1(jjm), yprimu1(jjm), rlatu2(jjm), yprimu2(jjm) 46 47 c 48 c ..... champs locaux ..... 49 c 50 51 REAL dzoom 52 REAL(KIND=8) ylat(jjp1), yprim(jjp1) 53 REAL(KIND=8) yuv 54 REAL(KIND=8) yt(0:nmax2) 55 REAL(KIND=8) fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2) 56 SAVE Ytprim, yt,Yf 57 REAL(KIND=8) Yf(0:nmax2),yypr(0:nmax2) 58 REAL(KIND=8) yvrai(jjp1), yprimm(jjp1),ylatt(jjp1) 59 REAL(KIND=8) pi,depi,pis2,epsilon,y0,pisjm 60 REAL(KIND=8) yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax 61 REAL(KIND=8) yfi,Yf1,ffdy 62 REAL(KIND=8) ypn,deply,y00 63 SAVE y00, deply 64 65 INTEGER i,j,it,ik,iter,jlat 66 INTEGER jpn,jjpn 67 SAVE jpn 68 REAL(KIND=8) a0,a1,a2,a3,yi2,heavyy0,heavyy0m 69 REAL(KIND=8) fa(0:nmax2),fb(0:nmax2) 70 REAL y0min,y0max 71 72 REAL(KIND=8) heavyside 73 74 pi = 2. * ASIN(1.) 75 depi = 2. * pi 76 pis2 = pi/2. 77 pisjm = pi/ REAL(jjm) 78 epsilon = 1.e-3 79 y0 = yzoomdeg * pi/180. 80 81 IF( dzooma.LT.1.) THEN 82 dzoom = dzooma * pi 83 ELSEIF( dzooma.LT. 12. ) THEN 84 WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug 85 ,menter et relancer ! ' 86 STOP 1 1 module fyhyp_m 2 3 IMPLICIT NONE 4 5 contains 6 7 SUBROUTINE fyhyp(rlatu, yyprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1) 8 9 ! From LMDZ4/libf/dyn3d/fyhyp.F, version 1.2, 2005/06/03 09:11:32 10 11 ! Author: P. Le Van, from analysis by R. Sadourny 12 13 ! Calcule les latitudes et dérivées dans la grille du GCM pour une 14 ! fonction f(y) à dérivée tangente hyperbolique. 15 16 ! Il vaut mieux avoir : grossismy * dzoom < pi / 2 17 18 use coefpoly_m, only: coefpoly 19 use nrtype, only: k8 20 21 include "dimensions.h" 22 ! for jjm 23 24 include "serre.h" 25 ! for clat, grossismy, dzoomy, tauy 26 27 REAL, intent(out):: rlatu(jjm + 1), yyprimu(jjm + 1) 28 REAL, intent(out):: rlatv(jjm) 29 real, intent(out):: rlatu2(jjm), yprimu2(jjm), rlatu1(jjm), yprimu1(jjm) 30 31 ! Local: 32 33 REAL(K8) champmin, champmax 34 INTEGER, PARAMETER:: nmax=30000, nmax2=2*nmax 35 REAL dzoom ! distance totale de la zone du zoom (en radians) 36 REAL(K8) ylat(jjm + 1), yprim(jjm + 1) 37 REAL(K8) yuv 38 REAL(K8), save:: yt(0:nmax2) 39 REAL(K8) fhyp(0:nmax2), beta 40 REAL(K8), save:: ytprim(0:nmax2) 41 REAL(K8) fxm(0:nmax2) 42 REAL(K8), save:: yf(0:nmax2) 43 REAL(K8) yypr(0:nmax2) 44 REAL(K8) yvrai(jjm + 1), yprimm(jjm + 1), ylatt(jjm + 1) 45 REAL(K8) pi, pis2, epsilon, y0, pisjm 46 REAL(K8) yo1, yi, ylon2, ymoy, yprimin 47 REAL(K8) yfi, yf1, ffdy 48 REAL(K8) ypn, deply, y00 49 SAVE y00, deply 50 51 INTEGER i, j, it, ik, iter, jlat 52 INTEGER jpn, jjpn 53 SAVE jpn 54 REAL(K8) a0, a1, a2, a3, yi2, heavyy0, heavyy0m 55 REAL(K8) fa(0:nmax2), fb(0:nmax2) 56 REAL y0min, y0max 57 58 REAL(K8) heavyside 59 60 !------------------------------------------------------------------- 61 62 print *, "Call sequence information: fyhyp" 63 64 pi = 2.*asin(1.) 65 pis2 = pi/2. 66 pisjm = pi/real(jjm) 67 epsilon = 1e-3 68 y0 = clat*pi/180. 69 dzoom = dzoomy*pi 70 print *, 'yzoom(rad), grossismy, tauy, dzoom (rad):' 71 print *, y0, grossismy, tauy, dzoom 72 73 DO i = 0, nmax2 74 yt(i) = -pis2 + real(i)*pi/nmax2 75 END DO 76 77 heavyy0m = heavyside(-y0) 78 heavyy0 = heavyside(y0) 79 y0min = 2.*y0*heavyy0m - pis2 80 y0max = 2.*y0*heavyy0 + pis2 81 82 fa = 999.999 83 fb = 999.999 84 85 DO i = 0, nmax2 86 IF (yt(i)<y0) THEN 87 fa(i) = tauy*(yt(i)-y0 + dzoom/2.) 88 fb(i) = (yt(i)-2.*y0*heavyy0m + pis2)*(y0-yt(i)) 89 ELSE IF (yt(i)>y0) THEN 90 fa(i) = tauy*(y0-yt(i) + dzoom/2.) 91 fb(i) = (2.*y0*heavyy0-yt(i) + pis2)*(yt(i)-y0) 92 END IF 93 94 IF (200.*fb(i)<-fa(i)) THEN 95 fhyp(i) = -1. 96 ELSE IF (200.*fb(i)<fa(i)) THEN 97 fhyp(i) = 1. 87 98 ELSE 88 dzoom = dzooma * pi/180. 89 ENDIF 90 91 WRITE(6,18) 92 WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)' 93 WRITE(6,24) y0,grossism,tau,dzoom 94 95 DO i = 0, nmax2 96 yt(i) = - pis2 + REAL(i)* pi /nmax2 97 ENDDO 98 99 heavyy0m = heavyside( -y0 ) 100 heavyy0 = heavyside( y0 ) 101 y0min = 2.*y0*heavyy0m - pis2 102 y0max = 2.*y0*heavyy0 + pis2 103 104 fa = 999.999 105 fb = 999.999 106 107 DO i = 0, nmax2 108 IF( yt(i).LT.y0 ) THEN 109 fa (i) = tau* (yt(i)-y0+dzoom/2. ) 110 fb(i) = (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) ) 111 ELSEIF ( yt(i).GT.y0 ) THEN 112 fa(i) = tau *(y0-yt(i)+dzoom/2. ) 113 fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 ) 114 ENDIF 115 116 IF( 200.* fb(i) .LT. - fa(i) ) THEN 117 fhyp ( i) = - 1. 118 ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN 119 fhyp ( i) = 1. 120 ELSE 121 fhyp(i) = TANH ( fa(i)/fb(i) ) 122 ENDIF 123 124 IF( yt(i).EQ.y0 ) fhyp(i) = 1. 125 IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1. 126 127 ENDDO 128 129 cc .... Calcul de beta .... 130 c 131 ffdy = 0. 132 133 DO i = 1, nmax2 134 ymoy = 0.5 * ( yt(i-1) + yt( i ) ) 135 IF( ymoy.LT.y0 ) THEN 136 fa(i)= tau * ( ymoy-y0+dzoom/2.) 137 fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy ) 138 ELSEIF ( ymoy.GT.y0 ) THEN 139 fa(i)= tau * ( y0-ymoy+dzoom/2. ) 140 fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 ) 141 ENDIF 142 143 IF( 200.* fb(i) .LT. - fa(i) ) THEN 144 fxm ( i) = - 1. 145 ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN 146 fxm ( i) = 1. 147 ELSE 148 fxm(i) = TANH ( fa(i)/fb(i) ) 149 ENDIF 150 IF( ymoy.EQ.y0 ) fxm(i) = 1. 151 IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1. 152 ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) ) 153 154 ENDDO 155 156 beta = ( grossism * ffdy - pi ) / ( ffdy - pi ) 157 158 IF( 2.*beta - grossism.LE. 0.) THEN 159 160 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 161 ,tine fyhyp est mauvaise ! ' 162 WRITE(6,*)'Modifier les valeurs de grossismy ,tauy ou dzoomy', 163 , ' et relancer ! *** ' 164 CALL ABORT_GCM("FYHYP", "", 1) 165 166 ENDIF 167 c 168 c ..... calcul de Ytprim ..... 169 c 170 171 DO i = 0, nmax2 172 Ytprim(i) = beta + ( grossism - beta ) * fhyp(i) 173 ENDDO 174 175 c ..... Calcul de Yf ........ 176 177 Yf(0) = - pis2 178 DO i = 1, nmax2 179 yypr(i) = beta + ( grossism - beta ) * fxm(i) 180 ENDDO 181 182 DO i=1,nmax2 183 Yf(i) = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) ) 184 ENDDO 185 186 c **************************************************************** 187 c 188 c ..... yuv = 0. si calcul des latitudes aux pts. U ..... 189 c ..... yuv = 0.5 si calcul des latitudes aux pts. V ..... 190 c 191 WRITE(6,18) 192 c 193 DO 5000 ik = 1,4 194 195 IF( ik.EQ.1 ) THEN 196 yuv = 0. 197 jlat = jjm + 1 198 ELSE IF ( ik.EQ.2 ) THEN 199 yuv = 0.5 200 jlat = jjm 201 ELSE IF ( ik.EQ.3 ) THEN 202 yuv = 0.25 203 jlat = jjm 204 ELSE IF ( ik.EQ.4 ) THEN 205 yuv = 0.75 206 jlat = jjm 207 ENDIF 208 c 209 yo1 = 0. 210 DO 1500 j = 1,jlat 211 yo1 = 0. 212 ylon2 = - pis2 + pisjm * ( REAL(j) + yuv -1.) 213 yfi = ylon2 214 c 215 DO 250 it = nmax2,0,-1 216 IF( yfi.GE.Yf(it)) GO TO 350 217 250 CONTINUE 218 it = 0 219 350 CONTINUE 220 221 yi = yt(it) 222 IF(it.EQ.nmax2) THEN 223 it = nmax2 -1 224 Yf(it+1) = pis2 225 ENDIF 226 c ................................................................. 227 c .... Interpolation entre yi(it) et yi(it+1) pour avoir Y(yi) 228 c ..... et Y'(yi) ..... 229 c ................................................................. 230 231 CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1), 232 , yt(it),yt(it+1) , a0,a1,a2,a3 ) 233 234 Yf1 = Yf(it) 235 Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi 236 237 DO 500 iter = 1,300 238 yi = yi - ( Yf1 - yfi )/ Yprimin 239 240 IF( ABS(yi-yo1).LE.epsilon) GO TO 550 241 yo1 = yi 242 yi2 = yi * yi 243 Yf1 = a0 + a1 * yi + a2 * yi2 + a3 * yi2 * yi 244 Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi2 245 500 CONTINUE 246 WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter 247 STOP 2 248 550 CONTINUE 249 c 250 Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi* yi 251 yprim(j) = pi / ( jjm * Yprimin ) 252 yvrai(j) = yi 253 254 1500 CONTINUE 255 256 DO j = 1, jlat -1 257 IF( yvrai(j+1). LT. yvrai(j) ) THEN 258 WRITE(6,*) ' PBS. avec rlat(',j+1,') plus petit que rlat(',j, 259 , ')' 260 STOP 3 261 ENDIF 262 ENDDO 263 264 WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2' 265 , ,' et pi/2 ' 266 c 267 IF( ik.EQ.1 ) THEN 268 ypn = pis2 269 DO j = jlat,1,-1 270 IF( yvrai(j).LE. ypn ) GO TO 1502 271 ENDDO 272 1502 CONTINUE 273 274 jpn = j 275 y00 = yvrai(jpn) 276 deply = pis2 - y00 277 ENDIF 278 279 DO j = 1, jjm +1 - jpn 280 ylatt (j) = -pis2 - y00 + yvrai(jpn+j-1) 281 yprimm(j) = yprim(jpn+j-1) 282 ENDDO 283 284 jjpn = jpn 285 IF( jlat.EQ. jjm ) jjpn = jpn -1 286 287 DO j = 1,jjpn 288 ylatt (j + jjm+1 -jpn) = yvrai(j) + deply 289 yprimm(j + jjm+1 -jpn) = yprim(j) 290 ENDDO 291 292 c *********** Fin de la reorganisation ************* 293 c 294 1600 CONTINUE 295 99 fhyp(i) = tanh(fa(i)/fb(i)) 100 END IF 101 102 IF (yt(i)==y0) fhyp(i) = 1. 103 IF (yt(i)==y0min .OR. yt(i)==y0max) fhyp(i) = -1. 104 END DO 105 106 ! Calcul de beta 107 108 ffdy = 0. 109 110 DO i = 1, nmax2 111 ymoy = 0.5*(yt(i-1) + yt(i)) 112 IF (ymoy<y0) THEN 113 fa(i) = tauy*(ymoy-y0 + dzoom/2.) 114 fb(i) = (ymoy-2.*y0*heavyy0m + pis2)*(y0-ymoy) 115 ELSE IF (ymoy>y0) THEN 116 fa(i) = tauy*(y0-ymoy + dzoom/2.) 117 fb(i) = (2.*y0*heavyy0-ymoy + pis2)*(ymoy-y0) 118 END IF 119 120 IF (200.*fb(i)<-fa(i)) THEN 121 fxm(i) = -1. 122 ELSE IF (200.*fb(i)<fa(i)) THEN 123 fxm(i) = 1. 124 ELSE 125 fxm(i) = tanh(fa(i)/fb(i)) 126 END IF 127 IF (ymoy==y0) fxm(i) = 1. 128 IF (ymoy==y0min .OR. yt(i)==y0max) fxm(i) = -1. 129 ffdy = ffdy + fxm(i)*(yt(i)-yt(i-1)) 130 END DO 131 132 beta = (grossismy*ffdy-pi)/(ffdy-pi) 133 134 IF (2. * beta - grossismy <= 0.) THEN 135 print *, 'Attention ! La valeur beta calculee dans la routine fyhyp ' & 136 // 'est mauvaise. Modifier les valeurs de grossismy, tauy ou ' & 137 // 'dzoomy et relancer.' 138 STOP 1 139 END IF 140 141 ! calcul de Ytprim 142 143 DO i = 0, nmax2 144 ytprim(i) = beta + (grossismy-beta)*fhyp(i) 145 END DO 146 147 ! Calcul de Yf 148 149 yf(0) = -pis2 150 DO i = 1, nmax2 151 yypr(i) = beta + (grossismy-beta)*fxm(i) 152 END DO 153 154 DO i = 1, nmax2 155 yf(i) = yf(i-1) + yypr(i)*(yt(i)-yt(i-1)) 156 END DO 157 158 ! yuv = 0. si calcul des latitudes aux pts. U 159 ! yuv = 0.5 si calcul des latitudes aux pts. V 160 161 loop_ik: DO ik = 1, 4 162 IF (ik==1) THEN 163 yuv = 0. 164 jlat = jjm + 1 165 ELSE IF (ik==2) THEN 166 yuv = 0.5 167 jlat = jjm 168 ELSE IF (ik==3) THEN 169 yuv = 0.25 170 jlat = jjm 171 ELSE IF (ik==4) THEN 172 yuv = 0.75 173 jlat = jjm 174 END IF 175 176 yo1 = 0. 296 177 DO j = 1, jlat 297 ylat(j) = ylatt( jlat +1 -j ) 298 yprim(j) = yprimm( jlat +1 -j ) 299 ENDDO 300 301 DO j = 1, jlat 302 yvrai(j) = ylat(j)*180./pi 303 ENDDO 304 305 IF( ik.EQ.1 ) THEN 306 c WRITE(6,18) 307 c WRITE(6,*) ' YLAT en U apres ( en deg. ) ' 308 c WRITE(6,68) (yvrai(j),j=1,jlat) 309 cc WRITE(6,*) ' YPRIM ' 310 cc WRITE(6,445) ( yprim(j),j=1,jlat) 311 312 DO j = 1, jlat 313 rrlatu(j) = ylat( j ) 314 yyprimu(j) = yprim( j ) 315 ENDDO 316 317 ELSE IF ( ik.EQ. 2 ) THEN 318 c WRITE(6,18) 319 c WRITE(6,*) ' YLAT en V apres ( en deg. ) ' 320 c WRITE(6,68) (yvrai(j),j=1,jlat) 321 cc WRITE(6,*)' YPRIM ' 322 cc WRITE(6,445) ( yprim(j),j=1,jlat) 323 324 DO j = 1, jlat 325 rrlatv(j) = ylat( j ) 326 yyprimv(j) = yprim( j ) 327 ENDDO 328 329 ELSE IF ( ik.EQ. 3 ) THEN 330 c WRITE(6,18) 331 c WRITE(6,*) ' YLAT en U + 0.75 apres ( en deg. ) ' 332 c WRITE(6,68) (yvrai(j),j=1,jlat) 333 cc WRITE(6,*) ' YPRIM ' 334 cc WRITE(6,445) ( yprim(j),j=1,jlat) 335 336 DO j = 1, jlat 337 rlatu2(j) = ylat( j ) 338 yprimu2(j) = yprim( j ) 339 ENDDO 340 341 ELSE IF ( ik.EQ. 4 ) THEN 342 c WRITE(6,18) 343 c WRITE(6,*) ' YLAT en U + 0.25 apres ( en deg. ) ' 344 c WRITE(6,68)(yvrai(j),j=1,jlat) 345 cc WRITE(6,*) ' YPRIM ' 346 cc WRITE(6,68) ( yprim(j),j=1,jlat) 347 348 DO j = 1, jlat 349 rlatu1(j) = ylat( j ) 350 yprimu1(j) = yprim( j ) 351 ENDDO 352 353 ENDIF 354 355 5000 CONTINUE 356 c 357 WRITE(6,18) 358 c 359 c ..... fin de la boucle do 5000 ..... 360 361 DO j = 1, jjm 362 ylat(j) = rrlatu(j) - rrlatu(j+1) 363 ENDDO 364 champmin = 1.e12 365 champmax = -1.e12 366 DO j = 1, jjm 367 champmin = MIN( champmin, ylat(j) ) 368 champmax = MAX( champmax, ylat(j) ) 369 ENDDO 370 champmin = champmin * 180./pi 371 champmax = champmax * 180./pi 372 373 24 FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3) 374 18 FORMAT(/) 375 68 FORMAT(1x,7f9.2) 376 377 RETURN 378 END 178 yo1 = 0. 179 ylon2 = -pis2 + pisjm*(real(j) + yuv-1.) 180 yfi = ylon2 181 182 it = nmax2 183 DO while (it >= 1 .and. yfi < yf(it)) 184 it = it - 1 185 END DO 186 187 yi = yt(it) 188 IF (it==nmax2) THEN 189 it = nmax2 - 1 190 yf(it + 1) = pis2 191 END IF 192 193 ! Interpolation entre yi(it) et yi(it + 1) pour avoir Y(yi) 194 ! et Y'(yi) 195 196 CALL coefpoly(yf(it), yf(it + 1), ytprim(it), ytprim(it + 1), & 197 yt(it), yt(it + 1), a0, a1, a2, a3) 198 199 yf1 = yf(it) 200 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi 201 202 iter = 1 203 DO 204 yi = yi - (yf1-yfi)/yprimin 205 IF (abs(yi-yo1)<=epsilon .or. iter == 300) exit 206 yo1 = yi 207 yi2 = yi*yi 208 yf1 = a0 + a1*yi + a2*yi2 + a3*yi2*yi 209 yprimin = a1 + 2.*a2*yi + 3.*a3*yi2 210 END DO 211 if (abs(yi-yo1) > epsilon) then 212 print *, 'Pas de solution.', j, ylon2 213 STOP 1 214 end if 215 216 yprimin = a1 + 2.*a2*yi + 3.*a3*yi*yi 217 yprim(j) = pi/(jjm*yprimin) 218 yvrai(j) = yi 219 END DO 220 221 DO j = 1, jlat - 1 222 IF (yvrai(j + 1)<yvrai(j)) THEN 223 print *, 'Problème avec rlat(', j + 1, ') plus petit que rlat(', & 224 j, ')' 225 STOP 1 226 END IF 227 END DO 228 229 print *, 'Reorganisation des latitudes pour avoir entre - pi/2 et pi/2' 230 231 IF (ik==1) THEN 232 ypn = pis2 233 DO j = jjm + 1, 1, -1 234 IF (yvrai(j)<=ypn) exit 235 END DO 236 237 jpn = j 238 y00 = yvrai(jpn) 239 deply = pis2 - y00 240 END IF 241 242 DO j = 1, jjm + 1 - jpn 243 ylatt(j) = -pis2 - y00 + yvrai(jpn + j-1) 244 yprimm(j) = yprim(jpn + j-1) 245 END DO 246 247 jjpn = jpn 248 IF (jlat==jjm) jjpn = jpn - 1 249 250 DO j = 1, jjpn 251 ylatt(j + jjm + 1-jpn) = yvrai(j) + deply 252 yprimm(j + jjm + 1-jpn) = yprim(j) 253 END DO 254 255 ! Fin de la reorganisation 256 257 DO j = 1, jlat 258 ylat(j) = ylatt(jlat + 1-j) 259 yprim(j) = yprimm(jlat + 1-j) 260 END DO 261 262 DO j = 1, jlat 263 yvrai(j) = ylat(j)*180./pi 264 END DO 265 266 IF (ik==1) THEN 267 DO j = 1, jjm + 1 268 rlatu(j) = ylat(j) 269 yyprimu(j) = yprim(j) 270 END DO 271 ELSE IF (ik==2) THEN 272 DO j = 1, jjm 273 rlatv(j) = ylat(j) 274 END DO 275 ELSE IF (ik==3) THEN 276 DO j = 1, jjm 277 rlatu2(j) = ylat(j) 278 yprimu2(j) = yprim(j) 279 END DO 280 ELSE IF (ik==4) THEN 281 DO j = 1, jjm 282 rlatu1(j) = ylat(j) 283 yprimu1(j) = yprim(j) 284 END DO 285 END IF 286 END DO loop_ik 287 288 DO j = 1, jjm 289 ylat(j) = rlatu(j) - rlatu(j + 1) 290 END DO 291 champmin = 1e12 292 champmax = -1e12 293 DO j = 1, jjm 294 champmin = min(champmin, ylat(j)) 295 champmax = max(champmax, ylat(j)) 296 END DO 297 champmin = champmin*180./pi 298 champmax = champmax*180./pi 299 300 DO j = 1, jjm 301 IF (rlatu1(j) <= rlatu2(j)) THEN 302 print *, 'Attention ! rlatu1 < rlatu2 ', rlatu1(j), rlatu2(j), j 303 STOP 13 304 ENDIF 305 306 IF (rlatu2(j) <= rlatu(j+1)) THEN 307 print *, 'Attention ! rlatu2 < rlatup1 ', rlatu2(j), rlatu(j+1), j 308 STOP 14 309 ENDIF 310 311 IF (rlatu(j) <= rlatu1(j)) THEN 312 print *, ' Attention ! rlatu < rlatu1 ', rlatu(j), rlatu1(j), j 313 STOP 15 314 ENDIF 315 316 IF (rlatv(j) <= rlatu2(j)) THEN 317 print *, ' Attention ! rlatv < rlatu2 ', rlatv(j), rlatu2(j), j 318 STOP 16 319 ENDIF 320 321 IF (rlatv(j) >= rlatu1(j)) THEN 322 print *, ' Attention ! rlatv > rlatu1 ', rlatv(j), rlatu1(j), j 323 STOP 17 324 ENDIF 325 326 IF (rlatv(j) >= rlatu(j)) THEN 327 print *, ' Attention ! rlatv > rlatu ', rlatv(j), rlatu(j), j 328 STOP 18 329 ENDIF 330 ENDDO 331 332 print *, 'Latitudes' 333 print 3, champmin, champmax 334 335 3 Format(1x, ' Au centre du zoom, la longueur de la maille est', & 336 ' d environ ', f0.2, ' degres ', /, & 337 ' alors que la maille en dehors de la zone du zoom est ', & 338 "d'environ ", f0.2, ' degres ') 339 340 END SUBROUTINE fyhyp 341 342 end module fyhyp_m -
trunk/LMDZ.COMMON/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90
r1422 r1441 31 31 INTEGER out_latudim,out_latvdim,out_dim(3) 32 32 INTEGER out_levdim 33 34 INTEGER, PARAMETER :: longcles = 2035 REAL clesphy0(longcles)36 33 37 34 INTEGER start(4),COUNT(4) … … 59 56 pa= 50000. 60 57 61 CALL conf_gcm( 99, .TRUE. , clesphy0)58 CALL conf_gcm( 99, .TRUE. ) 62 59 CALL iniconst 63 60 CALL inigeom -
trunk/LMDZ.COMMON/libf/dyn3d_common/inigeom.F
r1422 r1441 20 20 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, 21 21 . alphax,alphay,taux,tauy,transx,transy,pxo,pyo 22 use fxhyp_m, only: fxhyp 23 use fyhyp_m, only: fyhyp 22 24 23 25 IMPLICIT NONE … … 266 268 WRITE(6,*)'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' 267 269 268 CALL fxyhyper( clat, grossismy, dzoomy, tauy , 269 , clon, grossismx, dzoomx, taux , 270 , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 , 271 , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 ) 272 270 CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1) 271 CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) 273 272 274 273 ENDIF -
trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90
r1422 r1441 10 10 ! M.A Filiberti (04/2002) 11 11 ! 12 USE parallel_lmdz 13 USE Write_Field_p 14 USE Bands 12 USE parallel_lmdz, ONLY: ij_begin,ij_end,OMP_CHUNK,pole_nord,pole_sud,& 13 setdistrib 14 USE Write_Field_p, ONLY: WriteField_p 15 USE Bands, ONLY: jj_Nb_Caldyn,jj_Nb_vanleer 15 16 USE mod_hallo 16 17 USE Vampir -
trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F90
r1422 r1441 39 39 ! -metres du zoom avec celles lues sur le fichier start . 40 40 ! 41 LOGICAL etatinit42 INTEGER tapedef41 LOGICAL,INTENT(IN) :: etatinit 42 INTEGER,INTENT(IN) :: tapedef 43 43 44 44 ! Declarations : … … 48 48 include "comdissnew.h" 49 49 include "iniprint.h" 50 51 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique52 ! #include "clesphys.h"53 50 ! 54 51 ! … … 905 902 !Config Help = extension en longitude de la zone du zoom 906 903 !Config ( fraction de la zone totale) 907 dzoomx = 0. 0904 dzoomx = 0.2 908 905 CALL getin('dzoomx',dzoomx) 909 906 … … 913 910 !Config Help = extension en latitude de la zone du zoom 914 911 !Config ( fraction de la zone totale) 915 dzoomy = 0. 0912 dzoomy = 0.2 916 913 CALL getin('dzoomy',dzoomy) 917 914 -
trunk/LMDZ.COMMON/libf/dyn3dpar/covcont_p.F
r1019 r1441 1 1 SUBROUTINE covcont_p (klevel,ucov, vcov, ucont, vcont ) 2 USE parallel_lmdz 2 USE parallel_lmdz, ONLY: ij_begin,ij_end,OMP_CHUNK, 3 & pole_nord, pole_sud 3 4 IMPLICIT NONE 4 5 -
trunk/LMDZ.COMMON/libf/dyn3dpar/mod_hallo.F90
r1019 r1441 1 1 module mod_Hallo 2 USE parallel_lmdz 2 USE mod_const_mpi, ONLY: COMM_LMDZ,MPI_REAL_LMDZ 3 USE parallel_lmdz, ONLY: using_mpi, mpi_size, mpi_rank, omp_chunk, omp_rank, & 4 pole_nord, pole_sud, jj_begin, jj_end, & 5 jj_begin_para, jj_end_para 3 6 implicit none 4 7 logical,save :: use_mpi_alloc -
trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90
r1302 r1441 1 1 ! 2 ! $Id: parallel.F90 1 575 2011-09-21 13:57:48Z jghattas$2 ! $Id: parallel.F90 1810 2013-07-24 08:06:39Z emillour $ 3 3 ! 4 moduleparallel_lmdz4 MODULE parallel_lmdz 5 5 USE mod_const_mpi 6 6 #ifdef CPP_IOIPSL … … 31 31 integer, save :: omp_size 32 32 !$OMP THREADPRIVATE(omp_rank) 33 34 ! Ehouarn: add "dummy variables" (which are in dyn3d_mem/parallel_lmdz.F90) 35 ! so that calfis_loc compiles even if using dyn3dpar 36 integer,save :: jjb_u 37 integer,save :: jje_u 38 integer,save :: jjnb_u 39 integer,save :: jjb_v 40 integer,save :: jje_v 41 integer,save :: jjnb_v 42 43 integer,save :: ijb_u 44 integer,save :: ije_u 45 integer,save :: ijnb_u 46 47 integer,save :: ijb_v 48 integer,save :: ije_v 49 integer,save :: ijnb_v 33 50 34 51 contains … … 167 184 !Config Desc = taille des blocs openmp 168 185 !Config Def = 1 169 !Config Help = defini la taille des packets d'it �ration openmp186 !Config Help = defini la taille des packets d'it�ration openmp 170 187 !Config distribue a chaque tache lors de l'entree dans une 171 188 !Config boucle parallelisee … … 624 641 ! NewField(ij_be 625 642 626 end moduleparallel_lmdz643 end MODULE parallel_lmdz -
trunk/LMDZ.COMMON/libf/grid/dimension/makdim
r492 r1441 9 9 echo " $0 [im] [jm] lm" 10 10 echo " where im, jm and lm are the dimensions" 11 exit 11 exit 1 12 fi 13 14 if (($1 % 8 != 0)) && (( $# == 3 )) 15 then 16 echo "The number of longitudes must be a multiple of 8." 17 echo "See the files dyn3d/groupe.F and dyn3dpar/groupe_p.F." 18 exit 1 12 19 fi 13 20 -
trunk/LMDZ.COMMON/libf/misc/coefpoly_m.F90
r1440 r1441 1 ! 2 ! $Header$ 3 ! 4 SUBROUTINE coefpoly ( Xf1, Xf2, Xprim1, Xprim2, xtild1,xtild2 , 5 , a0,a1,a2,a3 ) 6 IMPLICIT NONE 7 c 8 c ... Auteur : P. Le Van ... 9 c 10 c 11 c Calcul des coefficients a0, a1, a2, a3 du polynome de degre 3 qui 12 c satisfait aux 4 equations suivantes : 1 module coefpoly_m 13 2 14 c a0 + a1*xtild1 + a2*xtild1*xtild1 + a3*xtild1*xtild1*xtild1 = Xf1 15 c a0 + a1*xtild2 + a2*xtild2*xtild2 + a3*xtild2*xtild2*xtild2 = Xf2 16 c a1 + 2.*a2*xtild1 + 3.*a3*xtild1*xtild1 = Xprim1 17 c a1 + 2.*a2*xtild2 + 3.*a3*xtild2*xtild2 = Xprim2 3 IMPLICIT NONE 18 4 19 c On en revient a resoudre un systeme de 4 equat.a 4 inconnues a0,a1,a2,a35 contains 20 6 21 REAL(KIND=8) Xf1, Xf2,Xprim1,Xprim2, xtild1,xtild2, xi 22 REAL(KIND=8) Xfout, Xprim 23 REAL(KIND=8) a1,a2,a3,a0, xtil1car, xtil2car,derr,x1x2car 7 SUBROUTINE coefpoly(xf1, xf2, xprim1, xprim2, xtild1, xtild2, a0, a1, a2, a3) 24 8 25 xtil1car = xtild1 * xtild1 26 xtil2car = xtild2 * xtild2 9 ! From LMDZ4/libf/dyn3d/coefpoly.F, version 1.1.1.1 2004/05/19 12:53:05 27 10 28 derr= 2. *(Xf2-Xf1)/( xtild1-xtild2)11 ! Author: P. Le Van 29 12 30 x1x2car = ( xtild1-xtild2)*(xtild1-xtild2) 13 ! Calcul des coefficients a0, a1, a2, a3 du polynôme de degré 3 qui 14 ! satisfait aux 4 équations suivantes : 31 15 32 a3 = (derr + Xprim1+Xprim2 )/x1x2car 33 a2 = ( Xprim1 - Xprim2 + 3.* a3 * ( xtil2car-xtil1car ) ) / 34 / ( 2.* ( xtild1 - xtild2 ) ) 16 ! a0 + a1 * xtild1 + a2 * xtild1**2 + a3 * xtild1**3 = Xf1 17 ! a0 + a1 * xtild2 + a2 * xtild2**2 + a3 * xtild2**3 = Xf2 18 ! a1 + 2. * a2 * xtild1 + 3. * a3 * xtild1**2 = Xprim1 19 ! a1 + 2. * a2 * xtild2 + 3. * a3 * xtild2**2 = Xprim2 35 20 36 a1 = Xprim1 -3.* a3 * xtil1car -2.* a2 * xtild137 a0 = Xf1 - a3 * xtild1* xtil1car -a2 * xtil1car - a1 *xtild121 ! (passe par les points (Xf(it), xtild(it)) et (Xf(it + 1), 22 ! xtild(it + 1)) 38 23 39 RETURN 40 END 24 ! On en revient à resoudre un système de 4 équations à 4 inconnues 25 ! a0, a1, a2, a3. 26 27 use nrtype, only: k8 28 29 REAL(K8), intent(in):: xf1, xf2, xprim1, xprim2, xtild1, xtild2 30 REAL(K8), intent(out):: a0, a1, a2, a3 31 32 ! Local: 33 REAL(K8) xtil1car, xtil2car, derr, x1x2car 34 35 !------------------------------------------------------------ 36 37 xtil1car = xtild1 * xtild1 38 xtil2car = xtild2 * xtild2 39 40 derr = 2. * (xf2-xf1)/(xtild1-xtild2) 41 42 x1x2car = (xtild1-xtild2) * (xtild1-xtild2) 43 44 a3 = (derr+xprim1+xprim2)/x1x2car 45 a2 = (xprim1-xprim2+3. * a3 * (xtil2car-xtil1car))/(2. * (xtild1-xtild2)) 46 47 a1 = xprim1 - 3. * a3 * xtil1car - 2. * a2 * xtild1 48 a0 = xf1 - a3 * xtild1 * xtil1car - a2 * xtil1car - a1 * xtild1 49 50 END SUBROUTINE coefpoly 51 52 end module coefpoly_m -
trunk/LMDZ.COMMON/libf/misc/wxios.F90
r1302 r1441 93 93 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 94 94 ! Routine d'initialisation !!!!!!!!!!!!! 95 ! A lancer juste apr ès mpi_init !!!!!!!!!!!!!95 ! A lancer juste après mpi_init !!!!!!!!!!!!! 96 96 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 97 97 … … 145 145 !Initialisation du contexte: 146 146 CALL xios_context_initialize(g_ctx_name, g_comm) 147 CALL xios_get_handle(g_ctx_name, xios_ctx) !R écupération147 CALL xios_get_handle(g_ctx_name, xios_ctx) !Récupération 148 148 CALL xios_set_current_context(xios_ctx) !Activation 149 149 g_ctx = xios_ctx … … 153 153 WRITE(lunout,*) " now call xios_solve_inheritance()" 154 154 ENDIF 155 !Une premi ère analyse des héritages:155 !Une première analyse des héritages: 156 156 CALL xios_solve_inheritance() 157 157 END SUBROUTINE wxios_context_init 158 158 159 159 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 160 ! Routine de param étrisation !!!!!!!!!!!!!!!!!!161 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 162 163 SUBROUTINE wxios_set_cal(pasdetemps, calendrier, annee, mois, jour, heure )164 IMPLICIT NONE 165 INCLUDE 'iniprint.h' 166 167 !Param ètres:160 ! Routine de paramétrisation !!!!!!!!!!!!!!!!!! 161 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 162 163 SUBROUTINE wxios_set_cal(pasdetemps, calendrier, annee, mois, jour, heure, ini_an, ini_mois, ini_jour, ini_heure) 164 IMPLICIT NONE 165 INCLUDE 'iniprint.h' 166 167 !Paramètres: 168 168 CHARACTER(len=*), INTENT(IN) :: calendrier 169 INTEGER, INTENT(IN) :: annee, mois, jour 170 REAL, INTENT(IN) :: pasdetemps, heure 169 INTEGER, INTENT(IN) :: annee, mois, jour, ini_an, ini_mois, ini_jour 170 REAL, INTENT(IN) :: pasdetemps, heure, ini_heure 171 171 172 172 !Variables: … … 181 181 mdtime = xios_time(0, 0, 0, 0, 0, pasdetemps) 182 182 183 !R églage du calendrier:183 !Réglage du calendrier: 184 184 SELECT CASE (calendrier) 185 185 CASE('earth_360d') … … 197 197 END SELECT 198 198 199 !Formatage de la date de départ: 200 WRITE(date, "(i4.4,'-',i2.2,'-',i2.2,' 00:00:00')") annee, mois, jour 201 202 IF (prt_level >= 10) WRITE(lunout,*) "wxios_set_cal: Initial time: ", date 203 204 CALL xios_set_context_attr_hdl(g_ctx, start_date= date) 199 !Formatage de la date d'origine: 200 WRITE(date, "(i4.4,'-',i2.2,'-',i2.2,' ',i2.2,':00:00')") annee, mois, jour, int(heure) 201 202 IF (prt_level >= 10) WRITE(lunout,*) "wxios_set_cal: Time origin: ", date 203 204 CALL xios_set_context_attr_hdl(g_ctx, time_origin = date) 205 206 !Formatage de la date de debut: 207 208 WRITE(date, "(i4.4,'-',i2.2,'-',i2.2,' ',i2.2,':00:00')") ini_an, ini_mois, ini_jour, int(ini_heure) 209 210 IF (prt_level >= 10) WRITE(lunout,*) "wxios_set_cal: Start date: ", date 211 212 CALL xios_set_context_attr_hdl(g_ctx, start_date = date) 205 213 206 214 !Et enfin,le pas de temps: … … 253 261 LOGICAL :: boool 254 262 255 !Masque pour les probl èmes de recouvrement MPI:263 !Masque pour les problèmes de recouvrement MPI: 256 264 LOGICAL :: mask(ni,nj) 257 265 258 !On r écupère le handle:266 !On récupère le handle: 259 267 CALL xios_get_domain_handle(dom_id, dom) 260 268 … … 285 293 286 294 CALL xios_is_defined_domain_attr_hdl(dom,ni_glo=boool) 287 !V érification:295 !Vérification: 288 296 IF (xios_is_valid_domain(dom_id)) THEN 289 297 IF (prt_level >= 10) WRITE(lunout,*) "wxios_domain_param: Domain initialized: ", trim(dom_id), boool … … 294 302 295 303 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 296 ! Pour d éclarer un axe vertical !!!!!!!!!!!!!!!304 ! Pour déclarer un axe vertical !!!!!!!!!!!!!!! 297 305 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 298 306 SUBROUTINE wxios_add_vaxis(axis_id, axis_size, axis_value) … … 315 323 ! axis_id=trim(axisgroup_id) 316 324 317 !On r écupère le groupe d'axes qui va bien:325 !On récupère le groupe d'axes qui va bien: 318 326 !CALL xios_get_axisgroup_handle(axisgroup_id, axgroup) 319 327 320 !On ajoute l'axe correspondant Ãce fichier:328 !On ajoute l'axe correspondant à ce fichier: 321 329 !CALL xios_add_axis(axgroup, ax, TRIM(ADJUSTL(axis_id))) 322 330 … … 327 335 CALL xios_set_axis_attr(trim(axis_id),size=axis_size,value=axis_value) 328 336 329 !V érification:337 !Vérification: 330 338 IF (xios_is_valid_axis(TRIM(ADJUSTL(axis_id)))) THEN 331 339 IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_vaxis: Axis created: ", TRIM(ADJUSTL(axis_id)) … … 338 346 339 347 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 340 ! Pour d éclarer un fichier !!!!!!!!!!!!!!!!!!!348 ! Pour déclarer un fichier !!!!!!!!!!!!!!!!!!! 341 349 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 342 350 SUBROUTINE wxios_add_file(fname, ffreq, flvl) … … 352 360 CHARACTER(len=100) :: nffreq 353 361 354 !On regarde si le fichier n'est pas d éfini par XML:362 !On regarde si le fichier n'est pas défini par XML: 355 363 IF (.NOT.xios_is_valid_file(fname)) THEN 356 !On cr ééle noeud:364 !On créé le noeud: 357 365 CALL xios_get_filegroup_handle("defile", x_fg) 358 366 CALL xios_add_file(x_fg, x_file, fname) 359 367 360 !On reformate la fr équence:368 !On reformate la fréquence: 361 369 CALL reformadate(ffreq, nffreq) 362 370 … … 376 384 ELSE 377 385 IF (prt_level >= 10) THEN 378 WRITE(lunout,*) "wxios_add_file: File ",trim(fname), " d éfined using XML."386 WRITE(lunout,*) "wxios_add_file: File ",trim(fname), " défined using XML." 379 387 ENDIF 380 388 ! Ehouarn: add an enable=.true. on top of xml definitions... why??? … … 384 392 385 393 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 386 ! Pour cr éer un champ !!!!!!!!!!!!!!!!!!!!394 ! Pour créer un champ !!!!!!!!!!!!!!!!!!!! 387 395 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 388 396 SUBROUTINE wxios_add_field(fieldname, fieldgroup, fieldlongname, fieldunit) … … 401 409 REAL(KIND=8) :: def 402 410 403 !La valeur par d éfaut des champs non définis:411 !La valeur par défaut des champs non définis: 404 412 def = nf90_fill_real 405 413 … … 414 422 !IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field: ",fieldname,fieldgroup, fieldlongname, fieldunit 415 423 416 !On rentre ses param ètres:424 !On rentre ses paramètres: 417 425 CALL xios_set_field_attr_hdl(field, standard_name=fieldlongname, unit=newunit, default_value=def) 418 426 IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field: Field ",trim(fieldname), "cree:" … … 422 430 423 431 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 424 ! Pour d éclarer un champ !!!!!!!!!!!!!!!!!425 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 426 SUBROUTINE wxios_add_field_to_file(fieldname, fdim, fid, fname, fieldlongname, fieldunit, field_level, op )432 ! Pour déclarer un champ !!!!!!!!!!!!!!!!! 433 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 434 SUBROUTINE wxios_add_field_to_file(fieldname, fdim, fid, fname, fieldlongname, fieldunit, field_level, op, nam_axvert) 427 435 IMPLICIT NONE 428 436 INCLUDE 'iniprint.h' … … 437 445 438 446 CHARACTER(len=20) :: axis_id ! Ehouarn: dangerous... 447 CHARACTER(len=20), INTENT(IN), OPTIONAL :: nam_axvert 439 448 CHARACTER(len=100) :: operation 440 449 TYPE(xios_file) :: f … … 451 460 axis_id="plev" 452 461 ENDIF 453 454 !on prépare le nom de l'opération: 462 463 IF (PRESENT(nam_axvert)) THEN 464 axis_id=nam_axvert 465 print*,'nam_axvert=',axis_id 466 ENDIF 467 468 !on prépare le nom de l'opération: 455 469 operation = reformaop(op) 456 470 … … 463 477 ENDIF 464 478 465 !On regarde si le champ à déjà été crééou non:479 !On regarde si le champ à déjà été créé ou non: 466 480 IF (xios_is_valid_field(fieldname) .AND. .NOT. g_field_name == fieldname) THEN 467 !Si ce champ existe via XML (ie, d ès le premier passage, ie g_field_name != fieldname) alors rien d'autre Ãfaire481 !Si ce champ existe via XML (ie, dès le premier passage, ie g_field_name != fieldname) alors rien d'autre à faire 468 482 IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field_to_file: Field ", trim(fieldname), "exists via XML" 469 483 g_flag_xml = .TRUE. … … 471 485 472 486 ELSE IF (.NOT. g_field_name == fieldname) THEN 473 !Si premier pssage et champ ind éfini, alors on le créé487 !Si premier pssage et champ indéfini, alors on le créé 474 488 475 489 IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field_to_file: Field ", trim(fieldname), "does not exist" 476 490 477 !On le cr éé:491 !On le créé: 478 492 CALL wxios_add_field(fieldname, fieldgroup, fieldlongname, fieldunit) 479 493 IF (xios_is_valid_field(fieldname)) THEN … … 487 501 488 502 IF (.NOT. g_flag_xml) THEN 489 !Champ existe d éjÃ, mais pas XML, alors on l'ajoute503 !Champ existe déjà, mais pas XML, alors on l'ajoute 490 504 !On ajoute le champ: 491 505 CALL xios_get_file_handle(fname, f) … … 497 511 498 512 499 !On rentre ses param ètres:513 !On rentre ses paramètres: 500 514 CALL xios_set_field_attr_hdl(field, level=field_level, enabled=.TRUE.) 501 515 … … 550 564 SUBROUTINE wxios_closedef() 551 565 CALL xios_close_context_definition() 552 CALL xios_update_calendar(0)566 ! CALL xios_update_calendar(0) 553 567 END SUBROUTINE wxios_closedef 554 568 … … 559 573 END MODULE wxios 560 574 #endif 561 -
trunk/LMDZ.COMMON/makelmdz_fcm
r1403 r1441 445 445 cd $LIBFGCM/grid/dimension 446 446 ./makdim $dim 447 if (($? != 0)) 448 then 449 exit 1 450 fi 451 447 452 cat $LIBFGCM/grid/dimensions.h 448 453 cd $LMDGCM
Note: See TracChangeset
for help on using the changeset viewer.
