Changeset 3876
- Timestamp:
- Aug 1, 2025, 9:37:44 PM (30 hours ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/dyn3d_common/disvert_noterre.F
r3874 r3876 3 3 4 4 c Auteur : F. Forget Y. Wanherdrick, P. Levan 5 c Nouvelle version 100% Mars !!6 c On l'utilise aussi pour Venus et Titan, legerment modifiee.7 5 8 6 #ifdef CPP_IOIPSL … … 24 22 c 25 23 c======================================================================= 26 c Discretisation verticale en coordonn�e hybride (ou sigma)27 c 24 c Vertical discretization for purely sigma or hybrid sigma-pressure 25 c coordinate depending on user provided flag "hybrid" 28 26 c======================================================================= 29 27 c … … 62 60 write(lunout,*) trim(modname),': hybrid=',hybrid 63 61 64 ! O uverture possible de fichiers typiquement E.T.62 ! Open file esasig.def, or if not there open z2sig.def 65 63 66 64 open(99,file="esasig.def",status='old',form='formatted', … … 73 71 74 72 c----------------------------------------------------------------------- 75 c cas 1 on lit les options dansesasig.def:73 c case 1 read from esasig.def: 76 74 c ---------------------------------------- 77 75 … … 128 126 129 127 c----------------------------------------------------------------------- 130 c cas 2 on lit les options dansz2sig.def:128 c case 2 read from z2sig.def: 131 129 c ---------------------------------------- 132 130 … … 147 145 read(99,*,iostat=ierr1) 148 146 if (ierr1 .eq. 0) then 149 write(lunout,*) 'ERROR: Expected nb of levels:', llm 150 call abort_gcm(modname, 151 & 'z2sig.def has more lines than expected', 1) 147 write(lunout,*) 'WARNING: More lines in z2sig.def' 148 & //' than expected from nb of levels:', llm 149 ! call abort_gcm(modname, 150 ! & 'z2sig.def has more lines than expected', 1) 152 151 endif 153 152 CLOSE(99) … … 178 177 179 178 c----------------------------------------------------------------------- 180 c .... C alculs de ap(l) et debp(l) ....179 c .... Compute ap(l) and bp(l) .... 181 180 c ......................................... 182 181 c 183 c ..... pa et preff sont lus sur les fichiers start par dynetat0 .....182 c ..... note that pa and preff are from start file read by dynetat0 184 183 c----------------------------------------------------------------------- 185 184 c … … 217 216 write(lunout,*) ap 218 217 219 c Calcul au milieu des couches : 220 c WARNING : le choix de placer le milieu des couches au niveau de 221 c pression interm�diaire est arbitraire et pourrait etre modifi�. 222 c Le calcul du niveau pour la derniere couche 223 c (on met la meme distance (en log pression) entre P(llm) 224 c et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est 225 c Specifique. Ce choix est sp�cifi� ici ET dans exner_milieu.F 218 c Compute value at the mid-layer : 219 c WARNING : there is an arbitrary decision here to put the 220 c middle of the layer half-way (pressure-wise) between boundaries 221 c In addition for the last layer (the top of which is at zero pressure) 222 c P(llm), we choose to put it with the same interval (in log(pressure)) 223 c from P(llm-1), than there is between between P(llm-1) and P(llm-2) 224 c This choice is quite specific and is done here AND in exner_milieu.F 225 c Both routines should definitely use the same assumptions. 226 226 227 227 DO l = 1, llm-1 … … 282 282 283 283 284 RETURN285 284 END 286 285 … … 292 291 c esasig.def/z2sig.def lors du passage en coordonnees hybrides 293 292 c F. Forget 2002 294 c Connaissant sig (niveaux "sigma" ou on veut mettre les couches)295 c L'objectif est de calculer newsig telle que293 c Knowing sig (the "sigma" levels where we want to put the layers) 294 c The goal is to compute newsig such that 296 295 c (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig 297 c Cela ne se r�soud pas analytiquement:298 c => on r�soud par iterration bourrine296 c Which cannot be done analyticaly: 297 c => need for an iterative procedure 299 298 c ---------------------------------------------- 300 c Information : whe re exp(1-1./x**2) become<< x299 c Information : when exp(1-1./x**2) becomes << x 301 300 c x exp(1-1./x**2) /x 302 301 c 1 1 … … 308 307 c 0.269 1.E-5 309 308 c 0.248 1.E-6 310 c => on peut utiliser newsig = sig*preff/pa sisig*preff/pa < 0.25309 c => one can thus use newsig = sig*preff/pa if sig*preff/pa < 0.25 311 310 312 311 … … 321 320 newsig= sig 322 321 else if (sig*preff/pa.ge.0.25) then 323 DO J=1,9999 ! nombre d''iteration max322 DO J=1,9999 ! overkill maximum number of iterations 324 323 F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig 325 324 c write(0,*) J, ' newsig =', newsig, ' F= ', F … … 331 330 newsig=(X2+newsig)*0.5 332 331 end if 333 c Test : on arete lorsque on approxime sig � moins de 0.01 m pr�s334 c (en pseudo altitude) :332 c Convergence is assumed if sig is within 0.01m (in pseudo-altitude) 333 c of the target value. 335 334 IF(abs(10.*log(F)).LT.1.E-5) goto 999 336 335 END DO … … 339 338 end if 340 339 999 continue 341 Return 340 342 341 END
Note: See TracChangeset
for help on using the changeset viewer.