Changeset 4368 for LMDZ6/branches/Ocean_skin/libf/dyn3d
- Timestamp:
- Dec 6, 2022, 12:01:16 AM (19 months ago)
- Location:
- LMDZ6/branches/Ocean_skin
- Files:
-
- 2 deleted
- 11 edited
- 3 copied
-
. (modified) (1 prop)
-
libf/dyn3d/advtrac.F90 (modified) (1 diff)
-
libf/dyn3d/check_isotopes.F (deleted)
-
libf/dyn3d/check_isotopes.F90 (copied) (copied from LMDZ6/trunk/libf/dyn3d/check_isotopes.F90)
-
libf/dyn3d/conf_gcm.F90 (modified) (10 diffs)
-
libf/dyn3d/dynetat0.F90 (copied) (copied from LMDZ6/trunk/libf/dyn3d/dynetat0.F90)
-
libf/dyn3d/dynetat0.f90 (deleted)
-
libf/dyn3d/dynredem.F90 (modified) (7 diffs)
-
libf/dyn3d/gcm.F90 (modified) (3 diffs)
-
libf/dyn3d/guide_mod.F90 (modified) (44 diffs)
-
libf/dyn3d/iniacademic.F90 (modified) (11 diffs)
-
libf/dyn3d/leapfrog.F (modified) (15 diffs)
-
libf/dyn3d/qminimum.F (modified) (11 diffs)
-
libf/dyn3d/replay3d.F90 (copied) (copied from LMDZ6/trunk/libf/dyn3d/replay3d.F90)
-
libf/dyn3d/vlsplt.F (modified) (26 diffs)
-
libf/dyn3d/vlspltqs.F (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Ocean_skin
- Property svn:mergeinfo changed
-
LMDZ6/branches/Ocean_skin/libf/dyn3d/advtrac.F90
r2622 r4368 1 1 ! $Id$ 2 2 3 SUBROUTINE advtrac(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk) 4 ! Auteur : F. Hourdin 5 ! 6 ! Modif. P. Le Van (20/12/97) 7 ! F. Codron (10/99) 8 ! D. Le Croller (07/2001) 9 ! M.A Filiberti (04/2002) 10 ! 11 USE infotrac, ONLY: nqtot, iadv,nqperes,ok_iso_verif 12 USE control_mod, ONLY: iapp_tracvl, day_step 13 USE comconst_mod, ONLY: dtvr 14 15 IMPLICIT NONE 16 ! 17 include "dimensions.h" 18 include "paramet.h" 19 include "comdissip.h" 20 include "comgeom2.h" 21 include "description.h" 22 include "iniprint.h" 23 24 !------------------------------------------------------------------- 25 ! Arguments 26 !------------------------------------------------------------------- 27 INTEGER,INTENT(OUT) :: iapptrac 28 REAL,INTENT(IN) :: pbaru(ip1jmp1,llm) 29 REAL,INTENT(IN) :: pbarv(ip1jm,llm) 30 REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) 31 REAL,INTENT(IN) :: masse(ip1jmp1,llm) 32 REAL,INTENT(IN) :: p( ip1jmp1,llmp1 ) 33 REAL,INTENT(IN) :: teta(ip1jmp1,llm) 34 REAL,INTENT(IN) :: pk(ip1jmp1,llm) 35 REAL,INTENT(OUT) :: flxw(ip1jmp1,llm) 36 !------------------------------------------------------------------- 37 ! Ajout PPM 38 !-------------------------------------------------------- 39 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) 40 !------------------------------------------------------------- 41 ! Variables locales 42 !------------------------------------------------------------- 43 44 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) 45 REAL massem(ip1jmp1,llm),zdp(ip1jmp1) 46 REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 47 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu 48 INTEGER iadvtr 49 INTEGER ij,l,iq,iiq 50 REAL zdpmin, zdpmax 51 EXTERNAL minmax 52 SAVE iadvtr, massem, pbaruc, pbarvc 53 DATA iadvtr/0/ 54 !---------------------------------------------------------- 55 ! Rajouts pour PPM 56 !---------------------------------------------------------- 57 INTEGER indice,n 58 REAL dtbon ! Pas de temps adaptatif pour que CFL<1 59 REAL CFLmaxz,aaa,bbb ! CFL maximum 60 REAL psppm(iim,jjp1) ! pression au sol 61 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) 62 REAL qppm(iim*jjp1,llm,nqtot) 63 REAL fluxwppm(iim,jjp1,llm) 64 REAL apppm(llmp1), bpppm(llmp1) 65 LOGICAL dum,fill 66 DATA fill/.true./ 67 DATA dum/.true./ 68 69 integer,save :: countcfl=0 70 real cflx(ip1jmp1,llm) 71 real cfly(ip1jm,llm) 72 real cflz(ip1jmp1,llm) 73 real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm) 74 75 IF(iadvtr.EQ.0) THEN 76 pbaruc(:,:)=0 77 pbarvc(:,:)=0 78 ENDIF 79 80 ! accumulation des flux de masse horizontaux 81 DO l=1,llm 82 DO ij = 1,ip1jmp1 83 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 84 ENDDO 85 DO ij = 1,ip1jm 86 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 87 ENDDO 88 ENDDO 89 90 ! selection de la masse instantannee des mailles avant le transport. 91 IF(iadvtr.EQ.0) THEN 92 3 #define DEBUG_IO 4 #undef DEBUG_IO 5 SUBROUTINE advtrac(pbaru, pbarv, p, masse,q,iapptrac,teta, flxw, pk) 6 ! Auteur : F. Hourdin 7 ! 8 ! Modif. P. Le Van (20/12/97) 9 ! F. Codron (10/99) 10 ! D. Le Croller (07/2001) 11 ! M.A Filiberti (04/2002) 12 ! 13 USE infotrac, ONLY: nqtot, tracers, isoCheck 14 USE control_mod, ONLY: iapp_tracvl, day_step 15 USE comconst_mod, ONLY: dtvr 16 17 IMPLICIT NONE 18 ! 19 include "dimensions.h" 20 include "paramet.h" 21 include "comdissip.h" 22 include "comgeom2.h" 23 include "description.h" 24 include "iniprint.h" 25 26 !--------------------------------------------------------------------------- 27 ! Arguments 28 !--------------------------------------------------------------------------- 29 INTEGER, INTENT(OUT) :: iapptrac 30 REAL, INTENT(IN) :: pbaru(ip1jmp1,llm) 31 REAL, INTENT(IN) :: pbarv(ip1jm, llm) 32 REAL, INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) 33 REAL, INTENT(IN) :: masse(ip1jmp1,llm) 34 REAL, INTENT(IN) :: p(ip1jmp1,llmp1 ) 35 REAL, INTENT(IN) :: teta(ip1jmp1,llm) 36 REAL, INTENT(IN) :: pk(ip1jmp1,llm) 37 REAL, INTENT(OUT) :: flxw(ip1jmp1,llm) 38 !--------------------------------------------------------------------------- 39 ! Ajout PPM 40 !--------------------------------------------------------------------------- 41 REAL :: massebx(ip1jmp1,llm), masseby(ip1jm,llm) 42 !--------------------------------------------------------------------------- 43 ! Variables locales 44 !--------------------------------------------------------------------------- 45 INTEGER :: ij, l, iq, iadv 46 ! REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu 47 REAL :: zdp(ip1jmp1), zdpmin, zdpmax 48 INTEGER, SAVE :: iadvtr=0 49 REAL, DIMENSION(ip1jmp1,llm) :: pbaruc, pbarug, massem, wg 50 REAL, DIMENSION(ip1jm, llm) :: pbarvc, pbarvg 51 EXTERNAL minmax 52 SAVE massem, pbaruc, pbarvc 53 !--------------------------------------------------------------------------- 54 ! Rajouts pour PPM 55 !--------------------------------------------------------------------------- 56 INTEGER indice, n 57 REAL :: dtbon ! Pas de temps adaptatif pour que CFL<1 58 REAL :: CFLmaxz, aaa, bbb ! CFL maximum 59 REAL, DIMENSION(iim,jjp1,llm) :: unatppm, vnatppm, fluxwppm 60 REAL :: qppm(iim*jjp1,llm,nqtot) 61 REAL :: psppm(iim,jjp1) ! pression au sol 62 REAL, DIMENSION(llmp1) :: apppm, bpppm 63 LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE. 64 65 INTEGER, SAVE :: countcfl=0 66 REAL, DIMENSION(ip1jmp1,llm) :: cflx, cflz 67 REAL, DIMENSION(ip1jm ,llm) :: cfly 68 REAL, DIMENSION(llm), SAVE :: cflxmax, cflymax, cflzmax 69 70 IF(iadvtr == 0) THEN 71 pbaruc(:,:)=0 72 pbarvc(:,:)=0 73 END IF 74 75 !--- Accumulation des flux de masse horizontaux 76 DO l=1,llm 77 DO ij = 1,ip1jmp1 78 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 79 END DO 80 DO ij = 1,ip1jm 81 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 82 END DO 83 END DO 84 85 !--- Selection de la masse instantannee des mailles avant le transport. 86 IF(iadvtr == 0) THEN 93 87 CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) 94 !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 95 ! 96 ENDIF 97 98 iadvtr = iadvtr+1 99 iapptrac = iadvtr 100 101 102 ! Test pour savoir si on advecte a ce pas de temps 103 IF ( iadvtr.EQ.iapp_tracvl ) THEN 104 105 !c .. Modif P.Le Van ( 20/12/97 ) .... 106 !c 107 108 ! traitement des flux de masse avant advection. 109 ! 1. calcul de w 110 ! 2. groupement des mailles pres du pole. 111 112 CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) 113 114 ! ... Flux de masse diaganostiques traceurs 115 flxw = wg / REAL(iapp_tracvl) 116 117 ! test sur l'eventuelle creation de valeurs negatives de la masse 118 DO l=1,llm-1 119 DO ij = iip2+1,ip1jm 120 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & 121 - pbarvg(ij-iip1,l) + pbarvg(ij,l) & 122 + wg(ij,l+1) - wg(ij,l) 123 ENDDO 124 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 125 DO ij = iip2,ip1jm 126 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 127 ENDDO 128 129 130 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 131 132 IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN 133 PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, & 134 ' MAX:', zdpmax 135 ENDIF 136 137 ENDDO 138 139 140 !------------------------------------------------------------------- 141 ! Calcul des criteres CFL en X, Y et Z 142 !------------------------------------------------------------------- 143 144 if (countcfl == 0. ) then 145 cflxmax(:)=0. 146 cflymax(:)=0. 147 cflzmax(:)=0. 148 endif 149 150 countcfl=countcfl+iapp_tracvl 151 cflx(:,:)=0. 152 cfly(:,:)=0. 153 cflz(:,:)=0. 154 do l=1,llm 155 do ij=iip2,ip1jm-1 156 if (pbarug(ij,l)>=0.) then 157 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) 158 else 159 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) 160 endif 161 enddo 162 enddo 163 do l=1,llm 164 do ij=iip2,ip1jm-1,iip1 165 cflx(ij+iip1,l)=cflx(ij,l) 166 enddo 167 enddo 168 169 do l=1,llm 170 do ij=1,ip1jm 171 if (pbarvg(ij,l)>=0.) then 172 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) 173 else 174 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) 175 endif 176 enddo 177 enddo 178 179 do l=2,llm 180 do ij=1,ip1jm 181 if (wg(ij,l)>=0.) then 182 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) 183 else 184 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) 185 endif 186 enddo 187 enddo 188 189 do l=1,llm 190 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) 191 cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) 192 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) 193 enddo 194 195 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 196 ! Par defaut, on sort le diagnostic des CFL tous les jours. 197 ! Si on veut le sortir a chaque pas d'advection en cas de plantage 198 ! if (countcfl==iapp_tracvl) then 199 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 200 if (countcfl==day_step) then 201 do l=1,llm 202 write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), & 203 cflzmax(l) 204 enddo 205 countcfl=0 206 endif 207 208 !------------------------------------------------------------------- 209 ! Advection proprement dite (Modification Le Croller (07/2001) 210 !------------------------------------------------------------------- 211 212 !---------------------------------------------------- 213 ! Calcul des moyennes basées sur la masse 214 !---------------------------------------------------- 215 call massbar(massem,massebx,masseby) 216 217 !----------------------------------------------------------- 218 ! Appel des sous programmes d'advection 219 !----------------------------------------------------------- 220 221 if (ok_iso_verif) then 222 write(*,*) 'advtrac 227' 223 call check_isotopes_seq(q,ip1jmp1,'advtrac 162') 224 endif !if (ok_iso_verif) then 225 226 do iq=1,nqperes 227 ! call clock(t_initial) 228 if(iadv(iq) == 0) cycle 229 ! ---------------------------------------------------------------- 230 ! Schema de Van Leer I MUSCL 231 ! ---------------------------------------------------------------- 232 if(iadv(iq).eq.10) THEN 233 ! CRisi: on fait passer tout q pour avoir acces aux fils 234 235 !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:) 236 call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) 237 238 ! ---------------------------------------------------------------- 239 ! Schema "pseudo amont" + test sur humidite specifique 240 ! pour la vapeur d'eau. F. Codron 241 ! ---------------------------------------------------------------- 242 else if(iadv(iq).eq.14) then 243 ! 244 !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:) 245 CALL vlspltqs( q, 2., massem, wg , & 246 pbarug,pbarvg,dtvr,p,pk,teta,iq) 247 248 ! ---------------------------------------------------------------- 249 ! Schema de Frederic Hourdin 250 ! ---------------------------------------------------------------- 251 else if(iadv(iq).eq.12) then 252 ! Pas de temps adaptatif 253 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 254 if (n.GT.1) then 255 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 256 dtvr,'n=',n 257 endif 258 do indice=1,n 259 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 260 end do 261 else if(iadv(iq).eq.13) then 262 ! Pas de temps adaptatif 263 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 264 if (n.GT.1) then 265 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 266 dtvr,'n=',n 267 endif 268 do indice=1,n 269 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 270 end do 271 ! ---------------------------------------------------------------- 272 ! Schema de pente SLOPES 273 ! ---------------------------------------------------------------- 274 else if (iadv(iq).eq.20) then 275 call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 276 277 ! ---------------------------------------------------------------- 278 ! Schema de Prather 279 ! ---------------------------------------------------------------- 280 else if (iadv(iq).eq.30) then 281 ! Pas de temps adaptatif 282 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 283 if (n.GT.1) then 284 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 285 dtvr,'n=',n 286 endif 287 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, & 288 n,dtbon) 289 290 ! ---------------------------------------------------------------- 291 ! Schemas PPM Lin et Rood 292 ! ---------------------------------------------------------------- 293 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. & 294 iadv(iq).LE.18)) then 295 296 ! Test sur le flux horizontal 297 ! Pas de temps adaptatif 298 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 299 if (n.GT.1) then 300 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 301 dtvr,'n=',n 302 endif 303 ! Test sur le flux vertical 304 CFLmaxz=0. 305 do l=2,llm 306 do ij=iip2,ip1jm 307 aaa=wg(ij,l)*dtvr/massem(ij,l) 308 CFLmaxz=max(CFLmaxz,aaa) 309 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 310 CFLmaxz=max(CFLmaxz,bbb) 311 enddo 312 enddo 313 if (CFLmaxz.GE.1) then 314 write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 315 endif 316 317 !----------------------------------------------------------- 318 ! Ss-prg interface LMDZ.4->PPM3d 319 !----------------------------------------------------------- 320 321 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & 322 apppm,bpppm,massebx,masseby,pbarug,pbarvg, & 323 unatppm,vnatppm,psppm) 324 325 do indice=1,n 326 !---------------------------------------------------------------- 327 ! VL (version PPM) horiz. et PPM vert. 328 !---------------------------------------------------------------- 329 if (iadv(iq).eq.11) then 330 ! Ss-prg PPM3d de Lin 331 call ppm3d(1,qppm(1,1,iq), & 332 psppm,psppm, & 333 unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, & 334 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 335 fill,dum,220.) 336 337 !------------------------------------------------------------- 338 ! Monotonic PPM 339 !------------------------------------------------------------- 340 else if (iadv(iq).eq.16) then 341 ! Ss-prg PPM3d de Lin 342 call ppm3d(1,qppm(1,1,iq), & 343 psppm,psppm, & 344 unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, & 345 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 346 fill,dum,220.) 347 !------------------------------------------------------------- 348 349 !------------------------------------------------------------- 350 ! Semi Monotonic PPM 351 !------------------------------------------------------------- 352 else if (iadv(iq).eq.17) then 353 ! Ss-prg PPM3d de Lin 354 call ppm3d(1,qppm(1,1,iq), & 355 psppm,psppm, & 356 unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, & 357 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 358 fill,dum,220.) 359 !------------------------------------------------------------- 360 361 !------------------------------------------------------------- 362 ! Positive Definite PPM 363 !------------------------------------------------------------- 364 else if (iadv(iq).eq.18) then 365 ! Ss-prg PPM3d de Lin 366 call ppm3d(1,qppm(1,1,iq), & 367 psppm,psppm, & 368 unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, & 369 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 370 fill,dum,220.) 371 !------------------------------------------------------------- 372 endif 373 enddo 374 !----------------------------------------------------------------- 375 ! Ss-prg interface PPM3d-LMDZ.4 376 !----------------------------------------------------------------- 377 call interpost(q(1,1,iq),qppm(1,1,iq)) 378 endif 379 !---------------------------------------------------------------------- 380 381 !----------------------------------------------------------------- 382 ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 383 ! et Nord j=1 384 !----------------------------------------------------------------- 385 386 ! call traceurpole(q(1,1,iq),massem) 387 388 ! calcul du temps cpu pour un schema donne 389 390 ! call clock(t_final) 391 !ym tps_cpu=t_final-t_initial 392 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu 393 394 end DO 395 396 if (ok_iso_verif) then 397 write(*,*) 'advtrac 402' 398 call check_isotopes_seq(q,ip1jmp1,'advtrac 397') 399 endif !if (ok_iso_verif) then 400 401 !------------------------------------------------------------------ 402 ! on reinitialise a zero les flux de masse cumules 403 !--------------------------------------------------- 404 iadvtr=0 405 406 ENDIF ! if iadvtr.EQ.iapp_tracvl 88 ! CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 89 END IF 90 91 iadvtr = iadvtr+1 92 iapptrac = iadvtr 93 94 !--- Test pour savoir si on advecte a ce pas de temps 95 IF(iadvtr /= iapp_tracvl) RETURN 96 97 ! .. Modif P.Le Van ( 20/12/97 ) .... 98 ! 99 ! traitement des flux de masse avant advection. 100 ! 1. calcul de w 101 ! 2. groupement des mailles pres du pole. 102 103 CALL groupe(massem, pbaruc, pbarvc, pbarug, pbarvg, wg) 104 105 !--- Flux de masse diaganostiques traceurs 106 flxw = wg / REAL(iapp_tracvl) 107 108 !--- Test sur l'eventuelle creation de valeurs negatives de la masse 109 DO l=1,llm-1 110 DO ij = iip2+1,ip1jm 111 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & 112 - pbarvg(ij-iip1,l) + pbarvg(ij,l) & 113 + wg(ij,l+1) - wg(ij,l) 114 END DO 115 ! ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ? 116 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 117 DO ij = iip2,ip1jm 118 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 119 END DO 120 121 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 122 123 IF(MAX(ABS(zdpmin),ABS(zdpmax)) > 0.5) & 124 WRITE(*,*)'WARNING DP/P l=',l,' MIN:',zdpmin,' MAX:', zdpmax 125 126 END DO 127 128 !------------------------------------------------------------------------- 129 ! Calcul des criteres CFL en X, Y et Z 130 !------------------------------------------------------------------------- 131 IF(countcfl == 0. ) then 132 cflxmax(:)=0. 133 cflymax(:)=0. 134 cflzmax(:)=0. 135 END IF 136 137 countcfl=countcfl+iapp_tracvl 138 cflx(:,:)=0. 139 cfly(:,:)=0. 140 cflz(:,:)=0. 141 DO l=1,llm 142 DO ij=iip2,ip1jm-1 143 IF(pbarug(ij,l)>=0.) then 144 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) 145 ELSE 146 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) 147 END IF 148 END DO 149 END DO 150 151 DO l=1,llm 152 DO ij=iip2,ip1jm-1,iip1 153 cflx(ij+iip1,l)=cflx(ij,l) 154 END DO 155 END DO 156 157 DO l=1,llm 158 DO ij=1,ip1jm 159 IF(pbarvg(ij,l)>=0.) then 160 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) 161 ELSE 162 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) 163 END IF 164 END DO 165 END DO 166 167 DO l=2,llm 168 DO ij=1,ip1jm 169 IF(wg(ij,l) >= 0.) THEN 170 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) 171 ELSE 172 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) 173 END IF 174 END DO 175 END DO 176 177 DO l=1,llm 178 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) 179 cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) 180 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) 181 END DO 182 183 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 184 ! Par defaut, on sort le diagnostic des CFL tous les jours. 185 ! Si on veut le sortir a chaque pas d'advection en cas de plantage 186 ! IF(countcfl==iapp_tracvl) then 187 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 188 IF(countcfl==day_step) then 189 DO l=1,llm 190 WRITE(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), cflzmax(l) 191 END DO 192 countcfl=0 193 END IF 194 195 !--------------------------------------------------------------------------- 196 ! Advection proprement dite (Modification Le Croller (07/2001) 197 !--------------------------------------------------------------------------- 198 199 !--------------------------------------------------------------------------- 200 ! Calcul des moyennes basees sur la masse 201 !--------------------------------------------------------------------------- 202 CALL massbar(massem,massebx,masseby) 203 204 #ifdef DEBUG_IO 205 CALL WriteField_u('massem',massem) 206 CALL WriteField_u('wg',wg) 207 CALL WriteField_u('pbarug',pbarug) 208 CALL WriteField_v('pbarvg',pbarvg) 209 CALL WriteField_u('p_tmp',p) 210 CALL WriteField_u('pk_tmp',pk) 211 CALL WriteField_u('teta_tmp',teta) 212 DO iq=1,nqtot 213 CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq)) 214 END DO 215 #endif 216 217 IF(isoCheck) WRITE(*,*) 'advtrac 227' 218 CALL check_isotopes_seq(q,ip1jmp1,'advtrac 162') 219 220 !------------------------------------------------------------------------- 221 ! Appel des sous programmes d'advection 222 !------------------------------------------------------------------------- 223 DO iq = 1, nqtot 224 ! CALL clock(t_initial) 225 IF(tracers(iq)%parent /= 'air') CYCLE 226 iadv = tracers(iq)%iadv 227 !----------------------------------------------------------------------- 228 SELECT CASE(iadv) 229 !----------------------------------------------------------------------- 230 CASE(0); CYCLE 231 !-------------------------------------------------------------------- 232 CASE(10) !--- Schema de Van Leer I MUSCL 233 !-------------------------------------------------------------------- 234 ! WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:) 235 CALL vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) 236 237 !-------------------------------------------------------------------- 238 CASE(14) !--- Schema "pseuDO amont" + test sur humidite specifique 239 !--- pour la vapeur d'eau. F. Codron 240 !-------------------------------------------------------------------- 241 ! WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:) 242 CALL vlspltqs(q,2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta,iq) 243 244 !-------------------------------------------------------------------- 245 CASE(12) !--- Schema de Frederic Hourdin 246 !-------------------------------------------------------------------- 247 CALL adaptdt(iadv,dtbon,n,pbarug,massem) ! pas de temps adaptatif 248 IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n 249 DO indice=1,n 250 CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 251 END DO 252 253 !-------------------------------------------------------------------- 254 CASE(13) !--- Pas de temps adaptatif 255 !-------------------------------------------------------------------- 256 CALL adaptdt(iadv,dtbon,n,pbarug,massem) 257 IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n 258 DO indice=1,n 259 CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 260 END DO 261 262 !-------------------------------------------------------------------- 263 CASE(20) !--- Schema de pente SLOPES 264 !-------------------------------------------------------------------- 265 CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 266 267 !-------------------------------------------------------------------- 268 CASE(30) !--- Schema de Prather 269 !-------------------------------------------------------------------- 270 ! Pas de temps adaptatif 271 CALL adaptdt(iadv,dtbon,n,pbarug,massem) 272 IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n 273 CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon) 274 275 !-------------------------------------------------------------------- 276 CASE(11,16,17,18) !--- Schemas PPM Lin et Rood 277 !-------------------------------------------------------------------- 278 ! Test sur le flux horizontal 279 CALL adaptdt(iadv,dtbon,n,pbarug,massem) ! pas de temps adaptatif 280 IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n 281 ! Test sur le flux vertical 282 CFLmaxz=0. 283 DO l=2,llm 284 DO ij=iip2,ip1jm 285 aaa=wg(ij,l)*dtvr/massem(ij,l) 286 CFLmaxz=max(CFLmaxz,aaa) 287 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 288 CFLmaxz=max(CFLmaxz,bbb) 289 END DO 290 END DO 291 IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 292 !---------------------------------------------------------------- 293 ! Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin) 294 !---------------------------------------------------------------- 295 CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & 296 apppm,bpppm,massebx,masseby,pbarug,pbarvg, & 297 unatppm,vnatppm,psppm) 298 299 !---------------------------------------------------------------- 300 DO indice=1,n !--- VL (version PPM) horiz. et PPM vert. 301 !---------------------------------------------------------------- 302 SELECT CASE(iadv) 303 !---------------------------------------------------------- 304 CASE(11) 305 !---------------------------------------------------------- 306 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & 307 2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) 308 !---------------------------------------------------------- 309 CASE(16) !--- Monotonic PPM 310 !---------------------------------------------------------- 311 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & 312 3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) 313 !---------------------------------------------------------- 314 CASE(17) !--- Semi monotonic PPM 315 !---------------------------------------------------------- 316 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & 317 4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.) 318 !---------------------------------------------------------- 319 CASE(18) !--- Positive Definite PPM 320 !---------------------------------------------------------- 321 CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, & 322 5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.) 323 END SELECT 324 !---------------------------------------------------------------- 325 END DO 326 !---------------------------------------------------------------- 327 ! Ss-prg interface PPM3d-LMDZ.4 328 !---------------------------------------------------------------- 329 CALL interpost(q(1,1,iq),qppm(1,1,iq)) 330 !---------------------------------------------------------------------- 331 END SELECT 332 !---------------------------------------------------------------------- 333 334 !---------------------------------------------------------------------- 335 ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1 336 !---------------------------------------------------------------------- 337 ! CALL traceurpole(q(1,1,iq),massem) 338 339 !--- Calcul du temps cpu pour un schema donne 340 ! CALL clock(t_final) 341 !ym tps_cpu=t_final-t_initial 342 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu 343 344 END DO 345 346 IF(isoCheck) WRITE(*,*) 'advtrac 402' 347 CALL check_isotopes_seq(q,ip1jmp1,'advtrac 397') 348 349 !------------------------------------------------------------------------- 350 ! on reinitialise a zero les flux de masse cumules 351 !------------------------------------------------------------------------- 352 iadvtr=0 407 353 408 354 END SUBROUTINE advtrac -
LMDZ6/branches/Ocean_skin/libf/dyn3d/conf_gcm.F90
r4013 r4368 116 116 CALL getin('calend', calend) 117 117 ! initialize year_len for aquaplanets and 1D 118 if (calend == 'earth_360d') then119 120 else if (calend == 'earth_365d') then121 122 else if (calend == 'earth_366d') then123 124 else125 126 endif118 IF (calend == 'earth_360d') THEN 119 year_len=360 120 ELSE IF (calend == 'earth_365d') THEN 121 year_len=365 122 ELSE IF (calend == 'earth_366d') THEN 123 year_len=366 124 ELSE 125 year_len=1 126 ENDIF 127 127 128 128 !Config Key = dayref … … 315 315 CALL getin('ngroup',ngroup) 316 316 317 318 317 ! top_bound sponge: only active if ok_strato=.true. and iflag_top_bound!=0 319 318 ! iflag_top_bound=0 for no sponge … … 322 321 iflag_top_bound=1 323 322 CALL getin('iflag_top_bound',iflag_top_bound) 323 IF (iflag_top_bound < 0 .or. iflag_top_bound > 2) & 324 call abort_gcm("conf_gcm", 'iflag_top_bound must be 0, 1 or 2', 1) 324 325 325 326 ! mode_top_bound : fields towards which sponge relaxation will be done: … … 394 395 ! ......... ( modif le 17/04/96 ) ......... 395 396 396 test_etatinit: IF (.not. etatinit) then397 test_etatinit: IF (.not. etatinit) THEN 397 398 !Config Key = clon 398 399 !Config Desc = centre du zoom, longitude … … 598 599 type_trac = 'lmdz' 599 600 CALL getin('type_trac',type_trac) 600 601 !Config Key = config_inca602 !Config Desc = Choix de configuration de INCA603 !Config Def = none604 !Config Help = Choix de configuration de INCA :605 !Config 'none' = sans INCA606 !Config 'chem' = INCA avec calcul de chemie607 !Config 'aero' = INCA avec calcul des aerosols608 config_inca = 'none'609 CALL getin('config_inca',config_inca)610 601 611 602 !Config Key = ok_dynzon … … 671 662 write(lunout,*)' offline = ', offline 672 663 write(lunout,*)' type_trac = ', type_trac 673 write(lunout,*)' config_inca = ', config_inca674 664 write(lunout,*)' ok_dynzon = ', ok_dynzon 675 665 write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 676 666 write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave 677 else667 ELSE 678 668 !Config Key = clon 679 669 !Config Desc = centre du zoom, longitude … … 795 785 CALL getin('type_trac',type_trac) 796 786 797 !Config Key = config_inca798 !Config Desc = Choix de configuration de INCA799 !Config Def = none800 !Config Help = Choix de configuration de INCA :801 !Config 'none' = sans INCA802 !Config 'chem' = INCA avec calcul de chemie803 !Config 'aero' = INCA avec calcul des aerosols804 config_inca = 'none'805 CALL getin('config_inca',config_inca)806 807 787 !Config Key = ok_dynzon 808 788 !Config Desc = sortie des transports zonaux dans la dynamique … … 875 855 876 856 write(lunout,*)' #########################################' 877 write(lunout,*)' Configuration des parametres de cel0' & 878 //'_limit: ' 857 write(lunout,*)' Configuration des parametres de cel0_limit: ' 879 858 write(lunout,*)' planet_type = ', planet_type 880 859 write(lunout,*)' calend = ', calend … … 912 891 write(lunout,*)' offline = ', offline 913 892 write(lunout,*)' type_trac = ', type_trac 914 write(lunout,*)' config_inca = ', config_inca915 893 write(lunout,*)' ok_dynzon = ', ok_dynzon 916 894 write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins … … 920 898 write(lunout,*)' ok_limit = ', ok_limit 921 899 write(lunout,*)' ok_etat0 = ', ok_etat0 900 write(lunout,*)' ok_guide = ', ok_guide 922 901 write(lunout,*)' read_orop = ', read_orop 923 endIF test_etatinit902 ENDIF test_etatinit 924 903 925 904 END SUBROUTINE conf_gcm -
LMDZ6/branches/Ocean_skin/libf/dyn3d/dynredem.F90
r4013 r4368 7 7 USE IOIPSL 8 8 #endif 9 USE infotrac 9 USE strings_mod, ONLY: maxlen 10 USE infotrac, ONLY: nqtot, tracers 10 11 USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL, & 11 12 NF90_CLOSE, NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER, & 12 13 NF90_64BIT_OFFSET 13 14 USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil 14 USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt,pa,preff, & 15 nivsig,nivsigs 15 USE comvert_mod, ONLY: ap, bp, presnivs, pa, preff, nivsig, nivsigs 16 16 USE comconst_mod, ONLY: cpp, daysec, dtvr, g, kappa, omeg, rad 17 17 USE logic_mod, ONLY: fxyhypb, ysinus … … 34 34 !=============================================================================== 35 35 ! Local variables: 36 INTEGER :: iq , l36 INTEGER :: iq 37 37 INTEGER, PARAMETER :: length=100 38 38 REAL :: tab_cntrl(length) !--- RUN PARAMETERS TABLE 39 39 ! For NetCDF: 40 CHARACTER(LEN= 30) :: unites40 CHARACTER(LEN=maxlen) :: unites 41 41 INTEGER :: indexID 42 42 INTEGER :: rlonuID, rlonvID, rlatuID, rlatvID 43 INTEGER :: sID, sigID, nID, vID,timID43 INTEGER :: sID, sigID, nID, timID 44 44 INTEGER :: yyears0, jjour0, mmois0 45 REAL :: z an0, zjulian, hours45 REAL :: zjulian, hours 46 46 !=============================================================================== 47 47 modname='dynredem0'; fil=fichnom … … 138 138 139 139 !--- Define fields saved later 140 WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')") ,&140 WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')") & 141 141 yyears0,mmois0,jjour0 142 142 CALL cre_var(nid,"temps","Temps de simulation",[timID],unites) … … 145 145 CALL cre_var(nid,"teta" ,"Temperature",[rlonvID,rlatuID,sID,timID]) 146 146 DO iq=1,nqtot 147 CALL cre_var(nid,t name(iq),ttext(iq),[rlonvID,rlatuID,sID,timID])147 CALL cre_var(nid,tracers(iq)%name,tracers(iq)%longName,[rlonvID,rlatuID,sID,timID]) 148 148 END DO 149 149 CALL cre_var(nid,"masse","Masse d air" ,[rlonvID,rlatuID,sID,timID]) … … 166 166 ! Purpose: Write the NetCDF restart file (append). 167 167 !------------------------------------------------------------------------------- 168 USE infotrac 168 USE strings_mod, ONLY: maxlen 169 USE infotrac, ONLY: nqtot, tracers, types_trac 169 170 USE control_mod 170 171 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_GET_VAR, NF90_INQ_VARID, & … … 192 193 !=============================================================================== 193 194 ! Local variables: 194 INTEGER :: l,iq, nid, vID, ierr, nid_trac, vID_trac195 INTEGER :: iq, nid, vID, ierr, nid_trac, vID_trac 195 196 INTEGER, SAVE :: nb=0 196 197 INTEGER, PARAMETER :: length=100 197 198 REAL :: tab_cntrl(length) ! tableau des parametres du run 198 CHARACTER(LEN= 256) :: var, dum199 CHARACTER(LEN=maxlen) :: var, dum 199 200 LOGICAL :: lread_inca 200 201 !=============================================================================== … … 227 228 !--- Tracers in file "start_trac.nc" (added by Anne) 228 229 lread_inca=.FALSE.; fil="start_trac.nc" 229 IF( type_trac=='inca' .OR. type_trac=='inco') INQUIRE(FILE=fil,EXIST=lread_inca)230 IF(ANY(types_trac=='inca') .OR. ANY(types_trac=='inco')) INQUIRE(FILE=fil,EXIST=lread_inca) 230 231 IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open") 231 232 232 233 !--- Save tracers 233 DO iq=1,nqtot; var= tname(iq); ierr=-1234 DO iq=1,nqtot; var=TRIM(tracers(iq)%name); ierr=-1 234 235 IF(lread_inca) THEN !--- Possibly read from "start_trac.nc" 235 236 fil="start_trac.nc" -
LMDZ6/branches/Ocean_skin/libf/dyn3d/gcm.F90
r3605 r4368 20 20 21 21 USE filtreg_mod 22 USE infotrac 22 USE infotrac, ONLY: nqtot, init_infotrac 23 23 USE control_mod 24 24 USE mod_const_mpi, ONLY: COMM_LMDZ … … 178 178 #ifdef CPP_IOIPSL 179 179 if (calend == 'earth_360d') then 180 call ioconf_calendar('360 d')180 call ioconf_calendar('360_day') 181 181 write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an' 182 182 else if (calend == 'earth_365d') then … … 205 205 ! Choix du nombre de traceurs et du schema pour l'advection 206 206 ! dans fichier traceur.def, par default ou via INCA 207 call in fotrac_init207 call init_infotrac 208 208 209 209 ! Allocation de la tableau q : champs advectes -
LMDZ6/branches/Ocean_skin/libf/dyn3d/guide_mod.F90
r4013 r4368 39 39 REAL, PRIVATE, SAVE :: lon_min_g,lon_max_g 40 40 REAL, PRIVATE, SAVE :: tau_lon,tau_lat 41 42 REAL, PRIVATE, SAVE :: plim_guide_BL 41 43 42 44 REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u,alpha_v … … 67 69 SUBROUTINE guide_init 68 70 71 use netcdf, only: nf90_noerr 69 72 USE control_mod, ONLY: day_step 70 73 USE serre_mod, ONLY: grossismx … … 113 116 CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie') 114 117 CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim') 118 CALL getpar('plim_guide_BL',85000.,plim_guide_BL,'BL top presnivs value') 119 115 120 116 121 ! Sauvegarde du for�age … … 169 174 if (ncidpl.eq.-99) then 170 175 rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl) 171 if (rcod.NE.NF _NOERR) THEN176 if (rcod.NE.NF90_NOERR) THEN 172 177 abort_message=' Nudging error -> no file apbp.nc' 173 178 CALL abort_gcm(modname,abort_message,1) … … 177 182 if (ncidpl.EQ.-99) then 178 183 rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl) 179 if (rcod.NE.NF _NOERR) THEN184 if (rcod.NE.NF90_NOERR) THEN 180 185 abort_message=' Nudging error -> no file P.nc' 181 186 CALL abort_gcm(modname,abort_message,1) … … 186 191 if (ncidpl.eq.-99) then 187 192 rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl) 188 if (rcod.NE.NF _NOERR) THEN193 if (rcod.NE.NF90_NOERR) THEN 189 194 CALL abort_gcm(modname, & 190 195 ' Nudging error -> no file u.nc',1) … … 195 200 if (ncidpl.eq.-99) then 196 201 rcod=nf90_open('v.nc',nf90_nowrite,ncidpl) 197 if (rcod.NE.NF _NOERR) THEN202 if (rcod.NE.NF90_NOERR) THEN 198 203 CALL abort_gcm(modname, & 199 204 ' Nudging error -> no file v.nc',1) … … 203 208 if (ncidpl.eq.-99) then 204 209 rcod=nf90_open('T.nc',nf90_nowrite,ncidpl) 205 if (rcod.NE.NF _NOERR) THEN210 if (rcod.NE.NF90_NOERR) THEN 206 211 CALL abort_gcm(modname, & 207 212 ' Nudging error -> no file T.nc',1) … … 211 216 if (ncidpl.eq.-99) then 212 217 rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl) 213 if (rcod.NE.NF _NOERR) THEN218 if (rcod.NE.NF90_NOERR) THEN 214 219 CALL abort_gcm(modname, & 215 220 ' Nudging error -> no file hur.nc',1) … … 220 225 endif 221 226 error=NF_INQ_DIMID(ncidpl,'LEVEL',rid) 222 IF (error.NE.NF _NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)223 IF (error.NE.NF _NOERR) THEN227 IF (error.NE.NF90_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid) 228 IF (error.NE.NF90_NOERR) THEN 224 229 CALL abort_gcm(modname,'Nudging: error reading pressure levels',1) 225 230 ENDIF … … 399 404 else 400 405 do l=1,llm 401 alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.406 alpha_pcor(l)=(1.+tanh(((plim_guide_BL-presnivs(l))/preff)/0.05))/2. 402 407 enddo 403 408 endif … … 1078 1083 SUBROUTINE guide_read(timestep) 1079 1084 1085 use netcdf, only: NF90_GET_VAR, nf90_noerr 1086 1080 1087 IMPLICIT NONE 1081 1088 1082 include "netcdf.inc"1083 1089 include "dimensions.h" 1084 1090 include "paramet.h" … … 1103 1109 if (first) then 1104 1110 ncidpl=-99 1105 write(*,*) ,trim(modname)//': opening nudging files '1111 write(*,*) trim(modname)//': opening nudging files ' 1106 1112 ! Niveaux de pression si non constants 1107 1113 if (guide_plevs.EQ.1) then 1108 write(*,*) ,trim(modname)//' Reading nudging on model levels'1114 write(*,*) trim(modname)//' Reading nudging on model levels' 1109 1115 rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl) 1110 IF (rcode.NE.NF _NOERR) THEN1116 IF (rcode.NE.NF90_NOERR) THEN 1111 1117 abort_message='Nudging: error -> no file apbp.nc' 1112 1118 CALL abort_gcm(modname,abort_message,1) 1113 1119 ENDIF 1114 1120 rcode = nf90_inq_varid(ncidpl, 'AP', varidap) 1115 IF (rcode.NE.NF _NOERR) THEN1121 IF (rcode.NE.NF90_NOERR) THEN 1116 1122 abort_message='Nudging: error -> no AP variable in file apbp.nc' 1117 1123 CALL abort_gcm(modname,abort_message,1) 1118 1124 ENDIF 1119 1125 rcode = nf90_inq_varid(ncidpl, 'BP', varidbp) 1120 IF (rcode.NE.NF _NOERR) THEN1126 IF (rcode.NE.NF90_NOERR) THEN 1121 1127 abort_message='Nudging: error -> no BP variable in file apbp.nc' 1122 1128 CALL abort_gcm(modname,abort_message,1) 1123 1129 ENDIF 1124 write(*,*) ,trim(modname)//' ncidpl,varidap',ncidpl,varidap1130 write(*,*) trim(modname)//' ncidpl,varidap',ncidpl,varidap 1125 1131 endif 1126 1132 … … 1128 1134 if (guide_plevs.EQ.2) then 1129 1135 rcode = nf90_open('P.nc', nf90_nowrite, ncidp) 1130 IF (rcode.NE.NF _NOERR) THEN1136 IF (rcode.NE.NF90_NOERR) THEN 1131 1137 abort_message='Nudging: error -> no file P.nc' 1132 1138 CALL abort_gcm(modname,abort_message,1) 1133 1139 ENDIF 1134 1140 rcode = nf90_inq_varid(ncidp, 'PRES', varidp) 1135 IF (rcode.NE.NF _NOERR) THEN1141 IF (rcode.NE.NF90_NOERR) THEN 1136 1142 abort_message='Nudging: error -> no PRES variable in file P.nc' 1137 1143 CALL abort_gcm(modname,abort_message,1) 1138 1144 ENDIF 1139 write(*,*) ,trim(modname)//' ncidp,varidp',ncidp,varidp1145 write(*,*) trim(modname)//' ncidp,varidp',ncidp,varidp 1140 1146 if (ncidpl.eq.-99) ncidpl=ncidp 1141 1147 endif … … 1144 1150 if (guide_u) then 1145 1151 rcode = nf90_open('u.nc', nf90_nowrite, ncidu) 1146 IF (rcode.NE.NF _NOERR) THEN1152 IF (rcode.NE.NF90_NOERR) THEN 1147 1153 abort_message='Nudging: error -> no file u.nc' 1148 1154 CALL abort_gcm(modname,abort_message,1) 1149 1155 ENDIF 1150 1156 rcode = nf90_inq_varid(ncidu, 'UWND', varidu) 1151 IF (rcode.NE.NF _NOERR) THEN1157 IF (rcode.NE.NF90_NOERR) THEN 1152 1158 abort_message='Nudging: error -> no UWND variable in file u.nc' 1153 1159 CALL abort_gcm(modname,abort_message,1) 1154 1160 ENDIF 1155 write(*,*) ,trim(modname)//' ncidu,varidu',ncidu,varidu1161 write(*,*) trim(modname)//' ncidu,varidu',ncidu,varidu 1156 1162 if (ncidpl.eq.-99) ncidpl=ncidu 1157 1163 … … 1175 1181 if (guide_v) then 1176 1182 rcode = nf90_open('v.nc', nf90_nowrite, ncidv) 1177 IF (rcode.NE.NF _NOERR) THEN1183 IF (rcode.NE.NF90_NOERR) THEN 1178 1184 abort_message='Nudging: error -> no file v.nc' 1179 1185 CALL abort_gcm(modname,abort_message,1) 1180 1186 ENDIF 1181 1187 rcode = nf90_inq_varid(ncidv, 'VWND', varidv) 1182 IF (rcode.NE.NF _NOERR) THEN1188 IF (rcode.NE.NF90_NOERR) THEN 1183 1189 abort_message='Nudging: error -> no VWND variable in file v.nc' 1184 1190 CALL abort_gcm(modname,abort_message,1) 1185 1191 ENDIF 1186 write(*,*) ,trim(modname)//' ncidv,varidv',ncidv,varidv1192 write(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv 1187 1193 if (ncidpl.eq.-99) ncidpl=ncidv 1188 1194 … … 1208 1214 if (guide_T) then 1209 1215 rcode = nf90_open('T.nc', nf90_nowrite, ncidt) 1210 IF (rcode.NE.NF _NOERR) THEN1216 IF (rcode.NE.NF90_NOERR) THEN 1211 1217 abort_message='Nudging: error -> no file T.nc' 1212 1218 CALL abort_gcm(modname,abort_message,1) 1213 1219 ENDIF 1214 1220 rcode = nf90_inq_varid(ncidt, 'AIR', varidt) 1215 IF (rcode.NE.NF _NOERR) THEN1221 IF (rcode.NE.NF90_NOERR) THEN 1216 1222 abort_message='Nudging: error -> no AIR variable in file T.nc' 1217 1223 CALL abort_gcm(modname,abort_message,1) 1218 1224 ENDIF 1219 write(*,*) ,trim(modname)//' ncidT,varidT',ncidt,varidt1225 write(*,*) trim(modname)//' ncidT,varidT',ncidt,varidt 1220 1226 if (ncidpl.eq.-99) ncidpl=ncidt 1221 1227 … … 1239 1245 if (guide_Q) then 1240 1246 rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ) 1241 IF (rcode.NE.NF _NOERR) THEN1247 IF (rcode.NE.NF90_NOERR) THEN 1242 1248 abort_message='Nudging: error -> no file hur.nc' 1243 1249 CALL abort_gcm(modname,abort_message,1) 1244 1250 ENDIF 1245 1251 rcode = nf90_inq_varid(ncidQ, 'RH', varidQ) 1246 IF (rcode.NE.NF _NOERR) THEN1252 IF (rcode.NE.NF90_NOERR) THEN 1247 1253 abort_message='Nudging: error -> no RH variable in file hur.nc' 1248 1254 CALL abort_gcm(modname,abort_message,1) 1249 1255 ENDIF 1250 write(*,*) ,trim(modname)//' ncidQ,varidQ',ncidQ,varidQ1256 write(*,*) trim(modname)//' ncidQ,varidQ',ncidQ,varidQ 1251 1257 if (ncidpl.eq.-99) ncidpl=ncidQ 1252 1258 … … 1270 1276 if ((guide_P).OR.(guide_modele)) then 1271 1277 rcode = nf90_open('ps.nc', nf90_nowrite, ncidps) 1272 IF (rcode.NE.NF _NOERR) THEN1278 IF (rcode.NE.NF90_NOERR) THEN 1273 1279 abort_message='Nudging: error -> no file ps.nc' 1274 1280 CALL abort_gcm(modname,abort_message,1) 1275 1281 ENDIF 1276 1282 rcode = nf90_inq_varid(ncidps, 'SP', varidps) 1277 IF (rcode.NE.NF _NOERR) THEN1283 IF (rcode.NE.NF90_NOERR) THEN 1278 1284 abort_message='Nudging: error -> no SP variable in file ps.nc' 1279 1285 CALL abort_gcm(modname,abort_message,1) 1280 1286 ENDIF 1281 write(*,*) ,trim(modname)//' ncidps,varidps',ncidps,varidps1287 write(*,*) trim(modname)//' ncidps,varidps',ncidps,varidps 1282 1288 endif 1283 1289 ! Coordonnee verticale … … 1285 1291 rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl) 1286 1292 IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl) 1287 write(*,*) ,trim(modname)//' ncidpl,varidpl',ncidpl,varidpl1293 write(*,*) trim(modname)//' ncidpl,varidpl',ncidpl,varidpl 1288 1294 endif 1289 1295 ! Coefs ap, bp pour calcul de la pression aux differents niveaux 1290 1296 if (guide_plevs.EQ.1) then 1291 #ifdef NC_DOUBLE 1292 status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc) 1293 status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc) 1294 #else 1295 status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc) 1296 status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc) 1297 #endif 1297 status=NF90_GET_VAR(ncidpl,varidap,apnc,[1],[nlevnc]) 1298 status=NF90_GET_VAR(ncidpl,varidbp,bpnc,[1],[nlevnc]) 1298 1299 ELSEIF (guide_plevs.EQ.0) THEN 1299 #ifdef NC_DOUBLE 1300 status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc) 1301 #else 1302 status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc) 1303 #endif 1300 status=NF90_GET_VAR(ncidpl,varidpl,apnc,[1],[nlevnc]) 1304 1301 !FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous 1305 1302 IF(convert_Pa) apnc=apnc*100.! conversion en Pascals … … 1326 1323 ! Pression 1327 1324 if (guide_plevs.EQ.2) then 1328 #ifdef NC_DOUBLE 1329 status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,pnat2) 1330 #else 1331 status=NF_GET_VARA_REAL(ncidp,varidp,start,count,pnat2) 1332 #endif 1325 status=NF90_GET_VAR(ncidp,varidp,pnat2,start,count) 1333 1326 IF (invert_y) THEN 1334 1327 ! PRINT*,"Invertion impossible actuellement" … … 1340 1333 ! Vent zonal 1341 1334 if (guide_u) then 1342 #ifdef NC_DOUBLE 1343 status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2) 1344 #else 1345 status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2) 1346 #endif 1335 status=NF90_GET_VAR(ncidu,varidu,unat2,start,count) 1347 1336 IF (invert_y) THEN 1348 1337 CALL invert_lat(iip1,jjp1,nlevnc,unat2) … … 1352 1341 ! Temperature 1353 1342 if (guide_T) then 1354 #ifdef NC_DOUBLE 1355 status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2) 1356 #else 1357 status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2) 1358 #endif 1343 status=NF90_GET_VAR(ncidt,varidt,tnat2,start,count) 1359 1344 IF (invert_y) THEN 1360 1345 CALL invert_lat(iip1,jjp1,nlevnc,tnat2) … … 1364 1349 ! Humidite 1365 1350 if (guide_Q) then 1366 #ifdef NC_DOUBLE 1367 status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2) 1368 #else 1369 status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2) 1370 #endif 1351 status=NF90_GET_VAR(ncidQ,varidQ,qnat2,start,count) 1371 1352 IF (invert_y) THEN 1372 1353 CALL invert_lat(iip1,jjp1,nlevnc,qnat2) … … 1378 1359 if (guide_v) then 1379 1360 count(2)=jjm 1380 #ifdef NC_DOUBLE 1381 status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2) 1382 #else 1383 status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2) 1384 #endif 1361 status=NF90_GET_VAR(ncidv,varidv,vnat2,start,count) 1385 1362 IF (invert_y) THEN 1386 1363 CALL invert_lat(iip1,jjm,nlevnc,vnat2) … … 1395 1372 count(3)=1 1396 1373 count(4)=0 1397 #ifdef NC_DOUBLE 1398 status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2) 1399 #else 1400 status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2) 1401 #endif 1374 status=NF90_GET_VAR(ncidps,varidps,psnat2,start,count) 1402 1375 IF (invert_y) THEN 1403 1376 CALL invert_lat(iip1,jjp1,1,psnat2) … … 1410 1383 SUBROUTINE guide_read2D(timestep) 1411 1384 1385 use netcdf, only: nf90_get_var, nf90_noerr 1386 1412 1387 IMPLICIT NONE 1413 1388 1414 include "netcdf.inc"1415 1389 include "dimensions.h" 1416 1390 include "paramet.h" … … 1443 1417 write(*,*)trim(modname)//' Reading nudging on model levels' 1444 1418 rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl) 1445 IF (rcode.NE.NF _NOERR) THEN1419 IF (rcode.NE.NF90_NOERR) THEN 1446 1420 abort_message='Nudging: error -> no file apbp.nc' 1447 1421 CALL abort_gcm(modname,abort_message,1) 1448 1422 ENDIF 1449 1423 rcode = nf90_inq_varid(ncidpl, 'AP', varidap) 1450 IF (rcode.NE.NF _NOERR) THEN1424 IF (rcode.NE.NF90_NOERR) THEN 1451 1425 abort_message='Nudging: error -> no AP variable in file apbp.nc' 1452 1426 CALL abort_gcm(modname,abort_message,1) 1453 1427 ENDIF 1454 1428 rcode = nf90_inq_varid(ncidpl, 'BP', varidbp) 1455 IF (rcode.NE.NF _NOERR) THEN1429 IF (rcode.NE.NF90_NOERR) THEN 1456 1430 abort_message='Nudging: error -> no BP variable in file apbp.nc' 1457 1431 CALL abort_gcm(modname,abort_message,1) … … 1462 1436 if (guide_plevs.EQ.2) then 1463 1437 rcode = nf90_open('P.nc', nf90_nowrite, ncidp) 1464 IF (rcode.NE.NF _NOERR) THEN1438 IF (rcode.NE.NF90_NOERR) THEN 1465 1439 abort_message='Nudging: error -> no file P.nc' 1466 1440 CALL abort_gcm(modname,abort_message,1) 1467 1441 ENDIF 1468 1442 rcode = nf90_inq_varid(ncidp, 'PRES', varidp) 1469 IF (rcode.NE.NF _NOERR) THEN1443 IF (rcode.NE.NF90_NOERR) THEN 1470 1444 abort_message='Nudging: error -> no PRES variable in file P.nc' 1471 1445 CALL abort_gcm(modname,abort_message,1) … … 1477 1451 if (guide_u) then 1478 1452 rcode = nf90_open('u.nc', nf90_nowrite, ncidu) 1479 IF (rcode.NE.NF _NOERR) THEN1453 IF (rcode.NE.NF90_NOERR) THEN 1480 1454 abort_message='Nudging: error -> no file u.nc' 1481 1455 CALL abort_gcm(modname,abort_message,1) 1482 1456 ENDIF 1483 1457 rcode = nf90_inq_varid(ncidu, 'UWND', varidu) 1484 IF (rcode.NE.NF _NOERR) THEN1458 IF (rcode.NE.NF90_NOERR) THEN 1485 1459 abort_message='Nudging: error -> no UWND variable in file u.nc' 1486 1460 CALL abort_gcm(modname,abort_message,1) … … 1492 1466 if (guide_v) then 1493 1467 rcode = nf90_open('v.nc', nf90_nowrite, ncidv) 1494 IF (rcode.NE.NF _NOERR) THEN1468 IF (rcode.NE.NF90_NOERR) THEN 1495 1469 abort_message='Nudging: error -> no file v.nc' 1496 1470 CALL abort_gcm(modname,abort_message,1) 1497 1471 ENDIF 1498 1472 rcode = nf90_inq_varid(ncidv, 'VWND', varidv) 1499 IF (rcode.NE.NF _NOERR) THEN1473 IF (rcode.NE.NF90_NOERR) THEN 1500 1474 abort_message='Nudging: error -> no VWND variable in file v.nc' 1501 1475 CALL abort_gcm(modname,abort_message,1) … … 1507 1481 if (guide_T) then 1508 1482 rcode = nf90_open('T.nc', nf90_nowrite, ncidt) 1509 IF (rcode.NE.NF _NOERR) THEN1483 IF (rcode.NE.NF90_NOERR) THEN 1510 1484 abort_message='Nudging: error -> no file T.nc' 1511 1485 CALL abort_gcm(modname,abort_message,1) 1512 1486 ENDIF 1513 1487 rcode = nf90_inq_varid(ncidt, 'AIR', varidt) 1514 IF (rcode.NE.NF _NOERR) THEN1488 IF (rcode.NE.NF90_NOERR) THEN 1515 1489 abort_message='Nudging: error -> no AIR variable in file T.nc' 1516 1490 CALL abort_gcm(modname,abort_message,1) … … 1522 1496 if (guide_Q) then 1523 1497 rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ) 1524 IF (rcode.NE.NF _NOERR) THEN1498 IF (rcode.NE.NF90_NOERR) THEN 1525 1499 abort_message='Nudging: error -> no file hur.nc' 1526 1500 CALL abort_gcm(modname,abort_message,1) 1527 1501 ENDIF 1528 1502 rcode = nf90_inq_varid(ncidQ, 'RH', varidQ) 1529 IF (rcode.NE.NF _NOERR) THEN1503 IF (rcode.NE.NF90_NOERR) THEN 1530 1504 abort_message='Nudging: error -> no RH,variable in file hur.nc' 1531 1505 CALL abort_gcm(modname,abort_message,1) … … 1537 1511 if ((guide_P).OR.(guide_modele)) then 1538 1512 rcode = nf90_open('ps.nc', nf90_nowrite, ncidps) 1539 IF (rcode.NE.NF _NOERR) THEN1513 IF (rcode.NE.NF90_NOERR) THEN 1540 1514 abort_message='Nudging: error -> no file ps.nc' 1541 1515 CALL abort_gcm(modname,abort_message,1) 1542 1516 ENDIF 1543 1517 rcode = nf90_inq_varid(ncidps, 'SP', varidps) 1544 IF (rcode.NE.NF _NOERR) THEN1518 IF (rcode.NE.NF90_NOERR) THEN 1545 1519 abort_message='Nudging: error -> no SP variable in file ps.nc' 1546 1520 CALL abort_gcm(modname,abort_message,1) … … 1556 1530 ! Coefs ap, bp pour calcul de la pression aux differents niveaux 1557 1531 if (guide_plevs.EQ.1) then 1558 #ifdef NC_DOUBLE 1559 status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc) 1560 status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc) 1561 #else 1562 status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc) 1563 status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc) 1564 #endif 1532 status=NF90_GET_VAR(ncidpl,varidap,apnc,[1],[nlevnc]) 1533 status=NF90_GET_VAR(ncidpl,varidbp,bpnc,[1],[nlevnc]) 1565 1534 elseif (guide_plevs.EQ.0) THEN 1566 #ifdef NC_DOUBLE 1567 status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc) 1568 #else 1569 status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc) 1570 #endif 1535 status=NF90_GET_VAR(ncidpl,varidpl,apnc,[1],[nlevnc]) 1571 1536 apnc=apnc*100.! conversion en Pascals 1572 1537 bpnc(:)=0. … … 1592 1557 ! Pression 1593 1558 if (guide_plevs.EQ.2) then 1594 #ifdef NC_DOUBLE 1595 status=NF_GET_VARA_DOUBLE(ncidp,varidp,start,count,zu) 1596 #else 1597 status=NF_GET_VARA_REAL(ncidp,varidp,start,count,zu) 1598 #endif 1559 status=NF90_GET_VAR(ncidp,varidp,zu,start,count) 1599 1560 DO i=1,iip1 1600 1561 pnat2(i,:,:)=zu(:,:) … … 1609 1570 ! Vent zonal 1610 1571 if (guide_u) then 1611 #ifdef NC_DOUBLE 1612 status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu) 1613 #else 1614 status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu) 1615 #endif 1572 status=NF90_GET_VAR(ncidu,varidu,zu,start,count) 1616 1573 DO i=1,iip1 1617 1574 unat2(i,:,:)=zu(:,:) … … 1626 1583 ! Temperature 1627 1584 if (guide_T) then 1628 #ifdef NC_DOUBLE 1629 status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu) 1630 #else 1631 status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu) 1632 #endif 1585 status=NF90_GET_VAR(ncidt,varidt,zu,start,count) 1633 1586 DO i=1,iip1 1634 1587 tnat2(i,:,:)=zu(:,:) … … 1643 1596 ! Humidite 1644 1597 if (guide_Q) then 1645 #ifdef NC_DOUBLE 1646 status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu) 1647 #else 1648 status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu) 1649 #endif 1598 status=NF90_GET_VAR(ncidQ,varidQ,zu,start,count) 1650 1599 DO i=1,iip1 1651 1600 qnat2(i,:,:)=zu(:,:) … … 1661 1610 if (guide_v) then 1662 1611 count(2)=jjm 1663 #ifdef NC_DOUBLE 1664 status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv) 1665 #else 1666 status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv) 1667 #endif 1612 status=NF90_GET_VAR(ncidv,varidv,zv,start,count) 1668 1613 DO i=1,iip1 1669 1614 vnat2(i,:,:)=zv(:,:) … … 1683 1628 count(3)=1 1684 1629 count(4)=0 1685 #ifdef NC_DOUBLE 1686 status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1)) 1687 #else 1688 status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1)) 1689 #endif 1630 status=NF90_GET_VAR(ncidps,varidps,zu(:,1),start,count) 1690 1631 DO i=1,iip1 1691 1632 psnat2(i,:)=zu(:,1) … … 1706 1647 USE comvert_mod, ONLY: presnivs 1707 1648 use netcdf95, only: nf95_def_var, nf95_put_var 1708 use netcdf, only: nf90_float 1649 use netcdf, only: nf90_float, nf90_def_var 1709 1650 1710 1651 IMPLICIT NONE … … 1748 1689 1749 1690 ! Creation des variables dimensions 1750 ierr=NF _DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)1751 ierr=NF _DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)1752 ierr=NF _DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)1753 ierr=NF _DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)1754 ierr=NF _DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)1755 ierr=NF _DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)1756 ierr=NF _DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)1757 ierr=NF _DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)1758 ierr=NF _DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)1691 ierr=NF90_DEF_VAR(nid,"LONU",NF90_FLOAT,id_lonu,vid_lonu) 1692 ierr=NF90_DEF_VAR(nid,"LONV",NF90_FLOAT,id_lonv,vid_lonv) 1693 ierr=NF90_DEF_VAR(nid,"LATU",NF90_FLOAT,id_latu,vid_latu) 1694 ierr=NF90_DEF_VAR(nid,"LATV",NF90_FLOAT,id_latv,vid_latv) 1695 ierr=NF90_DEF_VAR(nid,"LEVEL",NF90_FLOAT,id_lev,vid_lev) 1696 ierr=NF90_DEF_VAR(nid,"cu",NF90_FLOAT,(/id_lonu,id_latu/),vid_cu) 1697 ierr=NF90_DEF_VAR(nid,"cv",NF90_FLOAT,(/id_lonv,id_latv/),vid_cv) 1698 ierr=NF90_DEF_VAR(nid,"au",NF90_FLOAT,(/id_lonu,id_latu/),vid_au) 1699 ierr=NF90_DEF_VAR(nid,"av",NF90_FLOAT,(/id_lonv,id_latv/),vid_av) 1759 1700 call nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), & 1760 1701 varid_alpha_t) … … 1794 1735 ! Pressure (GCM) 1795 1736 dim4=(/id_lonv,id_latu,id_lev,id_tim/) 1796 ierr = NF _DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)1737 ierr = NF90_DEF_VAR(nid,"SP",NF90_FLOAT,dim4,varid) 1797 1738 ! Surface pressure (guidage) 1798 1739 IF (guide_P) THEN 1799 1740 dim3=(/id_lonv,id_latu,id_tim/) 1800 ierr = NF _DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)1741 ierr = NF90_DEF_VAR(nid,"ps",NF90_FLOAT,dim3,varid) 1801 1742 ENDIF 1802 1743 ! Zonal wind 1803 1744 IF (guide_u) THEN 1804 1745 dim4=(/id_lonu,id_latu,id_lev,id_tim/) 1805 ierr = NF _DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)1806 ierr = NF _DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)1807 ierr = NF _DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)1746 ierr = NF90_DEF_VAR(nid,"u",NF90_FLOAT,dim4,varid) 1747 ierr = NF90_DEF_VAR(nid,"ua",NF90_FLOAT,dim4,varid) 1748 ierr = NF90_DEF_VAR(nid,"ucov",NF90_FLOAT,dim4,varid) 1808 1749 ENDIF 1809 1750 ! Merid. wind 1810 1751 IF (guide_v) THEN 1811 1752 dim4=(/id_lonv,id_latv,id_lev,id_tim/) 1812 ierr = NF _DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)1813 ierr = NF _DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)1814 ierr = NF _DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)1753 ierr = NF90_DEF_VAR(nid,"v",NF90_FLOAT,dim4,varid) 1754 ierr = NF90_DEF_VAR(nid,"va",NF90_FLOAT,dim4,varid) 1755 ierr = NF90_DEF_VAR(nid,"vcov",NF90_FLOAT,dim4,varid) 1815 1756 ENDIF 1816 1757 ! Pot. Temperature 1817 1758 IF (guide_T) THEN 1818 1759 dim4=(/id_lonv,id_latu,id_lev,id_tim/) 1819 ierr = NF _DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)1760 ierr = NF90_DEF_VAR(nid,"teta",NF90_FLOAT,dim4,varid) 1820 1761 ENDIF 1821 1762 ! Specific Humidity 1822 1763 IF (guide_Q) THEN 1823 1764 dim4=(/id_lonv,id_latu,id_lev,id_tim/) 1824 ierr = NF _DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)1765 ierr = NF90_DEF_VAR(nid,"q",NF90_FLOAT,dim4,varid) 1825 1766 ENDIF 1826 1767 -
LMDZ6/branches/Ocean_skin/libf/dyn3d/iniacademic.F90
r4013 r4368 5 5 6 6 USE filtreg_mod, ONLY: inifilr 7 USE infotrac 7 USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName 8 8 USE control_mod, ONLY: day_step,planet_type 9 use exner_hyb_m, only: exner_hyb 10 use exner_milieu_m, only: exner_milieu 9 11 #ifdef CPP_IOIPSL 10 12 USE IOIPSL, ONLY: getin … … 14 16 #endif 15 17 USE Write_Field 16 use exner_hyb_m, only: exner_hyb17 use exner_milieu_m, only: exner_milieu18 18 USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm 19 19 USE logic_mod, ONLY: iflag_phys, read_start … … 21 21 USE temps_mod, ONLY: annee_ref, day_ini, day_ref 22 22 USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 23 USE readTracFiles_mod, ONLY: addPhase 23 24 24 25 ! Author: Frederic Hourdin original: 15/01/93 … … 61 62 real tetastrat ! potential temperature in the stratosphere, in K 62 63 real tetajl(jjp1,llm) 63 INTEGER i,j,l,lsup,ij 64 INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent 64 65 65 66 REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T … … 72 73 integer idum 73 74 74 REAL zdtvr 75 REAL zdtvr, tnat, alpha_ideal 75 76 76 77 character(len=*),parameter :: modname="iniacademic" 77 78 character(len=80) :: abort_message 78 79 79 80 80 ! Sanity check: verify that options selected by user are not incompatible … … 130 130 ps(:)=preff 131 131 endif 132 132 133 ! ground geopotential 133 134 phis(:)=0. 134 135 CALL pression ( ip1jmp1, ap, bp, ps, p ) 136 135 137 if (pressure_exner) then 136 138 CALL exner_hyb( ip1jmp1, ps, p, pks, pk) … … 193 195 CALL getin('tetanoise',tetanoise) 194 196 195 196 197 ! 2. Initialize fields towards which to relax 197 198 ! Friction … … 246 247 IF (.NOT. read_start) THEN 247 248 ! bulk initialization of temperature 248 249 249 IF (iflag_phys>10000) THEN 250 250 ! Particular case to impose a constant temperature T0=0.01*iflag_physx … … 253 253 teta(:,:)=tetarappel(:,:) 254 254 ENDIF 255 256 255 ! geopotential 257 256 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) … … 275 274 if (planet_type=="earth") then 276 275 ! Earth: first two tracers will be water 277 do i =1,nqtot278 if (i == 1) q(:,:,i)=1.e-10279 if (i == 2) q(:,:,i)=1.e-15280 if (i.gt.2) q(:,:,i)=0.276 do iq=1,nqtot 277 q(:,:,iq)=0. 278 IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10 279 IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15 281 280 282 281 ! CRisi: init des isotopes 283 282 ! distill de Rayleigh très simplifiée 284 if (ok_isotopes) then 285 if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then 286 q(:,:,i)=q(:,:,iqpere(i)) & 287 & *tnat(iso_num(i)) & 288 & *(q(:,:,iqpere(i))/30.e-3) & 289 & **(alpha_ideal(iso_num(i))-1) 290 endif 291 if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then 292 q(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i))) 293 endif 294 endif !if (ok_isotopes) then 295 283 iName = tracers(iq)%iso_iName 284 if (niso <= 0 .OR. iName <= 0) CYCLE 285 iPhase = tracers(iq)%iso_iPhase 286 iqParent = tracers(iq)%iqParent 287 IF(tracers(iq)%iso_iZone == 0) THEN 288 IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) & 289 CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1) 290 q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.) 291 ELSE 292 q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase)) 293 END IF 296 294 enddo 297 295 else … … 299 297 endif ! of if (planet_type=="earth") 300 298 301 if (ok_iso_verif) then 302 call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc') 303 endif !if (ok_iso_verif) then 299 call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc') 304 300 305 301 ! add random perturbation to temperature -
LMDZ6/branches/Ocean_skin/libf/dyn3d/leapfrog.F
r4013 r4368 11 11 use IOIPSL 12 12 #endif 13 USE infotrac, ONLY: nqtot, ok_iso_verif13 USE infotrac, ONLY: nqtot, isoCheck 14 14 USE guide_mod, ONLY : guide_main 15 15 USE write_field, ONLY: writefield … … 26 26 USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref, 27 27 & start_time,dt 28 USE strings_mod, ONLY: msg 28 29 29 30 IMPLICIT NONE … … 237 238 jH_cur = jH_cur - int(jH_cur) 238 239 239 if (ok_iso_verif) then 240 call check_isotopes_seq(q,ip1jmp1,'leapfrog 321') 241 endif !if (ok_iso_verif) then 240 call check_isotopes_seq(q,ip1jmp1,'leapfrog 321') 242 241 243 242 #ifdef CPP_IOIPSL … … 271 270 ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 ) 272 271 273 if (ok_iso_verif) then 274 call check_isotopes_seq(q,ip1jmp1,'leapfrog 400') 275 endif !if (ok_iso_verif) then 272 call check_isotopes_seq(q,ip1jmp1,'leapfrog 400') 276 273 277 274 2 CONTINUE ! Matsuno backward or leapfrog step begins here … … 324 321 325 322 326 if (ok_iso_verif) then 327 call check_isotopes_seq(q,ip1jmp1,'leapfrog 589') 328 endif !if (ok_iso_verif) then 323 call check_isotopes_seq(q,ip1jmp1,'leapfrog 589') 329 324 330 325 c----------------------------------------------------------------------- … … 345 340 c ------------------------------------------------------------- 346 341 347 if (ok_iso_verif) then 348 call check_isotopes_seq(q,ip1jmp1, 342 call check_isotopes_seq(q,ip1jmp1, 349 343 & 'leapfrog 686: avant caladvtrac') 350 endif !if (ok_iso_verif) then351 344 352 345 IF( forward. OR . leapf ) THEN … … 376 369 c ---------------------------------- 377 370 378 if (ok_iso_verif) then 379 write(*,*) 'leapfrog 720' 380 call check_isotopes_seq(q,ip1jmp1,'leapfrog 756') 381 endif !if (ok_iso_verif) then 371 CALL msg('720', modname, isoCheck) 372 call check_isotopes_seq(q,ip1jmp1,'leapfrog 756') 382 373 383 374 CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 , … … 385 376 ! $ finvmaold ) 386 377 387 if (ok_iso_verif) then 388 write(*,*) 'leapfrog 724' 389 call check_isotopes_seq(q,ip1jmp1,'leapfrog 762') 390 endif !if (ok_iso_verif) then 378 CALL msg('724', modname, isoCheck) 379 call check_isotopes_seq(q,ip1jmp1,'leapfrog 762') 391 380 392 381 c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd) … … 451 440 c+jld 452 441 453 c Diagnostique de conservation de l' énergie : initialisation442 c Diagnostique de conservation de l'energie : initialisation 454 443 IF (ip_ebil_dyn.ge.1 ) THEN 455 444 ztit='bil dyn' … … 498 487 499 488 c 500 c Diagnostique de conservation de l' énergie : difference489 c Diagnostique de conservation de l'energie : difference 501 490 IF (ip_ebil_dyn.ge.1 ) THEN 502 491 ztit='bil phys' … … 552 541 CALL massdair(p,masse) 553 542 554 if (ok_iso_verif) then 555 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196') 556 endif !if (ok_iso_verif) then 543 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196') 557 544 558 545 c----------------------------------------------------------------------- … … 639 626 c preparation du pas d'integration suivant ...... 640 627 641 if (ok_iso_verif) then 642 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509') 643 endif !if (ok_iso_verif) then 628 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509') 644 629 645 630 IF ( .NOT.purmats ) THEN … … 703 688 ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin)) 704 689 705 if (ok_iso_verif) then 706 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584') 707 endif !if (ok_iso_verif) then 690 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584') 708 691 709 692 c----------------------------------------------------------------------- … … 790 773 ELSE ! of IF (.not.purmats) 791 774 792 if (ok_iso_verif) then 793 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664') 794 endif !if (ok_iso_verif) then 775 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664') 795 776 796 777 c ........................................................ … … 817 798 ELSE ! of IF(forward) i.e. backward step 818 799 819 if (ok_iso_verif) then 820 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698') 821 endif !if (ok_iso_verif) then 800 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698') 822 801 823 802 IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN -
LMDZ6/branches/Ocean_skin/libf/dyn3d/qminimum.F
r2600 r4368 4 4 SUBROUTINE qminimum( q,nqtot,deltap ) 5 5 6 USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif 6 USE infotrac, ONLY: niso, ntiso,iqIsoPha, tracers 7 USE strings_mod, ONLY: strIdx 8 USE readTracFiles_mod, ONLY: addPhase 7 9 IMPLICIT none 8 10 c … … 16 18 REAL q(ip1jmp1,llm,nqtot), deltap(ip1jmp1,llm) 17 19 c 18 INTEGER iq_vap, iq_liq 19 PARAMETER ( iq_vap = 1 ) ! indice pour l'eau vapeur 20 PARAMETER ( iq_liq = 2 ) ! indice pour l'eau liquide 21 REAL seuil_vap, seuil_liq 22 PARAMETER ( seuil_vap = 1.0e-10 ) ! seuil pour l'eau vapeur 23 PARAMETER ( seuil_liq = 1.0e-11 ) ! seuil pour l'eau liquide 20 LOGICAL, SAVE :: first=.TRUE. 21 INTEGER, SAVE :: iq_vap, iq_liq ! indices pour l'eau vapeur/liquide 22 REAL, PARAMETER :: seuil_vap = 1.0e-10 ! seuil pour l'eau vapeur 23 REAL, PARAMETER :: seuil_liq = 1.0e-11 ! seuil pour l'eau liquide 24 24 c 25 25 c NB. ....( Il est souhaitable mais non obligatoire que les valeurs des … … 43 43 !INTEGER nb_pump 44 44 INTEGER ixt 45 46 IF(first) THEN 47 iq_vap = strIdx(tracers(:)%name, addPhase('H2O', 'g')) 48 iq_liq = strIdx(tracers(:)%name, addPhase('H2O', 'l')) 49 first = .FALSE. 50 END IF 45 51 c 46 52 c Quand l'eau liquide est trop petite (ou negative), on prend … … 49 55 c 50 56 51 if (ok_iso_verif) then 52 call check_isotopes_seq(q,ip1jmp1,'qminimum 52') 53 endif !if (ok_iso_verif) then 57 call check_isotopes_seq(q,ip1jmp1,'qminimum 52') 54 58 55 59 zx_defau_diag(:,:,:)=0.0 … … 59 63 if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then 60 64 61 if (ok_isotopes) then 62 zx_defau_diag(i,k,iq_liq)=AMAX1 65 if (niso > 0) zx_defau_diag(i,k,iq_liq)=AMAX1 63 66 : ( seuil_liq - q(i,k,iq_liq), 0.0 ) 64 endif !if (ok_isotopes) then65 67 66 68 q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq … … 80 82 if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then 81 83 82 if (ok_isotopes) then 83 zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 ) 84 endif !if (ok_isotopes) then 84 if (niso > 0) 85 & zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 ) 85 86 86 87 q(i,k-1,iq) = q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) * … … 110 111 111 112 !write(*,*) 'qminimum 128' 112 if ( ok_isotopes) then113 if (niso > 0) then 113 114 ! CRisi: traiter de même les traceurs d'eau 114 115 ! Mais il faut les prendre à l'envers pour essayer de conserver la … … 130 131 if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 131 132 ! on ajoute la vapeur en k 132 do ixt=1,nt raciso133 q(i,k,iq iso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))134 : 135 : *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)133 do ixt=1,ntiso 134 q(i,k,iqIsoPha(ixt,iq_vap))=q(i,k,iqIsoPha(ixt,iq_vap)) 135 : +zx_defau_diag(i,k,iq_vap) 136 : *q(i,k-1,iqIsoPha(ixt,iq_vap))/q_follow(i,k-1,iq_vap) 136 137 137 138 ! et on la retranche en k-1 138 q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap)) 139 q(i,k-1,iqIsoPha(ixt,iq_vap))= 140 : q(i,k-1,iqIsoPha(ixt,iq_vap)) 139 141 : -zx_defau_diag(i,k,iq_vap) 140 142 : *deltap(i,k)/deltap(i,k-1) 141 : *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap) 143 : *q(i,k-1,iqIsoPha(ixt,iq_vap)) 144 : /q_follow(i,k-1,iq_vap) 142 145 143 146 enddo !do ixt=1,niso … … 151 154 enddo !do k=2,llm 152 155 153 if (ok_iso_verif) then 154 call check_isotopes_seq(q,ip1jmp1,'qminimum 168') 155 endif !if (ok_iso_verif) then 156 call check_isotopes_seq(q,ip1jmp1,'qminimum 168') 156 157 157 158 … … 163 164 164 165 ! on ajoute eau liquide en k en k 165 do ixt=1,nt raciso166 q(i,k,iq iso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq))166 do ixt=1,ntiso 167 q(i,k,iqIsoPha(ixt,iq_liq))=q(i,k,iqIsoPha(ixt,iq_liq)) 167 168 : +zx_defau_diag(i,k,iq_liq) 168 : *q(i,k,iq iso(ixt,iq_vap))/q_follow(i,k,iq_vap)169 : *q(i,k,iqIsoPha(ixt,iq_vap))/q_follow(i,k,iq_vap) 169 170 ! et on la retranche à la vapeur en k 170 q(i,k,iq iso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))171 q(i,k,iqIsoPha(ixt,iq_vap))=q(i,k,iqIsoPha(ixt,iq_vap)) 171 172 : -zx_defau_diag(i,k,iq_liq) 172 : *q(i,k,iq iso(ixt,iq_vap))/q_follow(i,k,iq_vap)173 : *q(i,k,iqIsoPha(ixt,iq_vap))/q_follow(i,k,iq_vap) 173 174 enddo !do ixt=1,niso 174 175 q_follow(i,k,iq_liq)= q_follow(i,k,iq_liq) … … 180 181 enddo !do k=2,llm 181 182 182 if (ok_iso_verif) then 183 call check_isotopes_seq(q,ip1jmp1,'qminimum 197') 184 endif !if (ok_iso_verif) then 183 call check_isotopes_seq(q,ip1jmp1,'qminimum 197') 185 184 186 endif !if ( ok_isotopes) then185 endif !if (niso > 0) then 187 186 !write(*,*) 'qminimum 188' 188 187 -
LMDZ6/branches/Ocean_skin/libf/dyn3d/vlsplt.F
r4013 r4368 4 4 5 5 SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq) 6 USE infotrac, ONLY: nqtot, nqdesc,iqfils6 USE infotrac, ONLY: nqtot,tracers 7 7 c 8 8 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 38 38 c --------- 39 39 c 40 INTEGER i,ij,l,j,ii 41 INTEGER ijlqmin,iqmin,jqmin,lqmin 42 c 43 REAL zm(ip1jmp1,llm,nqtot),newmasse 40 INTEGER ij,l 41 c 42 REAL zm(ip1jmp1,llm,nqtot) 44 43 REAL mu(ip1jmp1,llm) 45 44 REAL mv(ip1jm,llm) 46 45 REAL mw(ip1jmp1,llm+1) 47 REAL zq(ip1jmp1,llm,nqtot),zz 48 REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm) 49 REAL second,temps0,temps1,temps2,temps3 50 REAL ztemps1,ztemps2,ztemps3 46 REAL zq(ip1jmp1,llm,nqtot) 51 47 REAL zzpbar, zzw 52 LOGICAL testcpu53 SAVE testcpu54 SAVE temps1,temps2,temps355 INTEGER iminn,imaxx56 48 INTEGER ifils,iq2 ! CRisi 57 49 58 50 REAL qmin,qmax 59 51 DATA qmin,qmax/0.,1.e33/ 60 DATA testcpu/.false./61 DATA temps1,temps2,temps3/0.,0.,0./62 63 52 64 53 zzpbar = 0.5 * pdt … … 83 72 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1) 84 73 85 if (nqdesc(iq).gt.0) then 86 do ifils=1,nqdesc(iq) 87 iq2=iqfils(ifils,iq) 88 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 89 enddo 90 endif !if (nqfils(iq).gt.0) then 74 do ifils=1,tracers(iq)%nqDescen 75 iq2=tracers(iq)%iqDescen(ifils) 76 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 77 enddo 91 78 92 79 cprint*,'Entree vlx1' … … 122 109 ENDDO 123 110 ! CRisi: aussi pour les fils 124 if (nqdesc(iq).gt.0) then 125 do ifils=1,nqdesc(iq) 126 iq2=iqfils(ifils,iq) 111 do ifils=1,tracers(iq)%nqDescen 112 iq2=tracers(iq)%iqDescen(ifils) 127 113 DO l=1,llm 128 DO ij=1,ip1jmp1129 q(ij,l,iq2)=zq(ij,l,iq2)130 ENDDO131 DO ij=1,ip1jm+1,iip1114 DO ij=1,ip1jmp1 115 q(ij,l,iq2)=zq(ij,l,iq2) 116 ENDDO 117 DO ij=1,ip1jm+1,iip1 132 118 q(ij+iim,l,iq2)=q(ij,l,iq2) 133 ENDDO119 ENDDO 134 120 ENDDO 135 enddo !do ifils=1,nqdesc(iq) 136 endif ! if (nqdesc(iq).gt.0) then 121 enddo 137 122 138 123 RETURN 139 124 END 140 125 RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq) 141 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi142 & qperemin,masseqmin,ratiomin! MVals et CRisi126 USE infotrac, ONLY : nqtot,tracers, ! CRisi 127 & min_qParent,min_qMass,min_ratio ! MVals et CRisi 143 128 144 129 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 161 146 c ---------- 162 147 REAL masse(ip1jmp1,llm,nqtot),pente_max 163 REAL u_m( ip1jmp1,llm ) ,pbarv( iip1,jjm,llm)148 REAL u_m( ip1jmp1,llm ) 164 149 REAL q(ip1jmp1,llm,nqtot) 165 REAL w(ip1jmp1,llm)166 150 INTEGER iq ! CRisi 167 151 c … … 173 157 c 174 158 REAL new_m,zu_m,zdum(ip1jmp1,llm) 175 REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1) 159 c REAL sigu(ip1jmp1) 160 REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1) 176 161 REAL zz(ip1jmp1) 177 162 REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm) … … 182 167 INTEGER ifils,iq2 ! CRisi 183 168 184 Logical extremum,first,testcpu 185 SAVE first,testcpu 186 187 REAL SSUM 188 REAL temps0,temps1,temps2,temps3,temps4,temps5,second 189 SAVE temps0,temps1,temps2,temps3,temps4,temps5 190 191 REAL z1,z2,z3 192 193 DATA first,testcpu/.true.,.false./ 194 195 IF(first) THEN 196 temps1=0. 197 temps2=0. 198 temps3=0. 199 temps4=0. 200 temps5=0. 201 first=.false. 202 ENDIF 169 Logical first 170 SAVE first 171 DATA first/.true./ 203 172 204 173 c calcul de la pente a droite et a gauche de la maille … … 436 405 ENDDO 437 406 ENDIF ! n0.gt.0 438 9999 continue407 c9999 continue 439 408 440 409 … … 449 418 ! CRisi: appel récursif de l'advection sur les fils. 450 419 ! Il faut faire ça avant d'avoir mis à jour q et masse 451 !write(*,*) 'vlsplt 326: iq,nq fils(iq)=',iq,nqfils(iq)420 !write(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen 452 421 453 if (nqdesc(iq).gt.0) then 454 do ifils=1,nqdesc(iq) 455 iq2=iqfils(ifils,iq) 456 DO l=1,llm 422 do ifils=1,tracers(iq)%nqDescen 423 iq2=tracers(iq)%iqDescen(ifils) 424 DO l=1,llm 457 425 DO ij=iip2,ip1jm 458 ! On a besoin de q et masse seulement entre iip2 et ip1jm459 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)460 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)461 !Mvals: veiller a ce qu'on n'ait pas de denominateur nul462 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)463 if (q(ij,l,iq).gt.qperemin) then464 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)465 else466 Ratio(ij,l,iq2)=ratiomin467 endif426 ! On a besoin de q et masse seulement entre iip2 et ip1jm 427 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 428 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 429 !Mvals: veiller a ce qu'on n'ait pas de denominateur nul 430 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 431 if (q(ij,l,iq).gt.min_qParent) then 432 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 433 else 434 Ratio(ij,l,iq2)=min_ratio 435 endif 468 436 enddo 469 enddo 470 enddo !do ifils=1,nqdesc(iq) 471 do ifils=1,nqfils(iq) 472 iq2=iqfils(ifils,iq) 473 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 474 enddo !do ifils=1,nqfils(iq) 475 endif !if (nqfils(iq).gt.0) then 437 enddo 438 enddo 439 do ifils=1,tracers(iq)%nqChildren 440 iq2=tracers(iq)%iqDescen(ifils) 441 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 442 enddo 476 443 ! end CRisi 477 444 … … 482 449 DO ij=iip2+1,ip1jm 483 450 !MVals: veiller a ce qu'on ait pas de denominateur nul 484 new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),m asseqmin)451 new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),min_qMass) 485 452 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ 486 453 & u_mq(ij-1,l)-u_mq(ij,l)) … … 498 465 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 499 466 ! puis on boucle en longitude 500 if (nqfils(iq).gt.0) then 501 do ifils=1,nqdesc(iq) 502 iq2=iqfils(ifils,iq) 503 DO l=1,llm 467 do ifils=1,tracers(iq)%nqDescen 468 iq2=tracers(iq)%iqDescen(ifils) 469 DO l=1,llm 504 470 DO ij=iip2+1,ip1jm 505 471 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) … … 508 474 q(ij-iim,l,iq2)=q(ij,l,iq2) 509 475 enddo ! DO ij=ijb+iip1-1,ije,iip1 510 enddo !DO l=1,llm 511 enddo !do ifils=1,nqdesc(iq) 512 endif !if (nqfils(iq).gt.0) then 476 enddo !DO l=1,llm 477 enddo 513 478 514 479 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) … … 519 484 END 520 485 RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq) 521 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi522 & qperemin,masseqmin,ratiomin! MVals et CRisi486 USE infotrac, ONLY : nqtot,tracers, ! CRisi 487 & min_qParent,min_qMass,min_ratio ! MVals et CRisi 523 488 c 524 489 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 544 509 REAL masse(ip1jmp1,llm,nqtot),pente_max 545 510 REAL masse_adv_v( ip1jm,llm) 546 REAL q(ip1jmp1,llm,nqtot) , dq( ip1jmp1,llm)511 REAL q(ip1jmp1,llm,nqtot) 547 512 INTEGER iq ! CRisi 548 513 c … … 553 518 c 554 519 REAL airej2,airejjm,airescb(iim),airesch(iim) 555 REAL dyq(ip1jmp1,llm),dyqv(ip1jm) ,zdvm(ip1jmp1,llm)520 REAL dyq(ip1jmp1,llm),dyqv(ip1jm) 556 521 REAL adyqv(ip1jm),dyqmax(ip1jmp1) 557 522 REAL qbyv(ip1jm,llm) 558 523 559 REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs 524 REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs 525 c REAL appn apps 560 526 c REAL newq,oldmasse 561 Logical extremum,first,testcpu 562 REAL temps0,temps1,temps2,temps3,temps4,temps5,second 563 SAVE temps0,temps1,temps2,temps3,temps4,temps5 564 SAVE first,testcpu 527 LOGICAL first 528 SAVE first 565 529 566 530 REAL convpn,convps,convmpn,convmps … … 578 542 REAL SSUM 579 543 580 DATA first,testcpu/.true.,.false./ 581 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./ 544 DATA first/.true./ 582 545 583 546 !write(*,*) 'vly 578: entree, iq=',iq … … 778 741 ! CRisi: appel récursif de l'advection sur les fils. 779 742 ! Il faut faire ça avant d'avoir mis à jour q et masse 780 !write(*,*) 'vly 689: iq,nq fils(iq)=',iq,nqfils(iq)743 !write(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen 781 744 782 if (nqfils(iq).gt.0) then 783 do ifils=1,nqdesc(iq) 784 iq2=iqfils(ifils,iq) 785 DO l=1,llm 786 DO ij=1,ip1jmp1 787 ! attention, chaque fils doit avoir son masseq, sinon, le 1er 788 ! fils ecrase le masseq de ses freres. 789 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 790 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 791 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 792 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin) 793 if (q(ij,l,iq).gt.qperemin) then 794 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 795 else 796 Ratio(ij,l,iq2)=ratiomin 797 endif 745 do ifils=1,tracers(iq)%nqDescen 746 iq2=tracers(iq)%iqDescen(ifils) 747 DO l=1,llm 748 DO ij=1,ip1jmp1 749 ! attention, chaque fils doit avoir son masseq, sinon, le 1er 750 ! fils ecrase le masseq de ses freres. 751 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 752 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 753 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 754 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 755 if (q(ij,l,iq).gt.min_qParent) then 756 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 757 else 758 Ratio(ij,l,iq2)=min_ratio 759 endif 798 760 enddo 799 enddo 800 enddo !do ifils=1,nqdesc(iq) 801 802 do ifils=1,nqfils(iq) 803 iq2=iqfils(ifils,iq) 804 call vly(Ratio,pente_max,masseq,qbyv,iq2) 805 enddo !do ifils=1,nqfils(iq) 806 endif !if (nqfils(iq).gt.0) then 761 enddo 762 enddo 763 764 do ifils=1,tracers(iq)%nqDescen 765 iq2=tracers(iq)%iqDescen(ifils) 766 call vly(Ratio,pente_max,masseq,qbyv,iq2) 767 enddo 807 768 808 769 DO l=1,llm … … 872 833 873 834 ! retablir les fils en rapport de melange par rapport a l'air: 874 if (nqfils(iq).gt.0) then 875 do ifils=1,nqdesc(iq) 876 iq2=iqfils(ifils,iq) 877 DO l=1,llm 835 do ifils=1,tracers(iq)%nqDescen 836 iq2=tracers(iq)%iqDescen(ifils) 837 DO l=1,llm 878 838 DO ij=1,ip1jmp1 879 839 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 880 840 enddo 881 enddo 882 enddo !do ifils=1,nqdesc(iq) 883 endif !if (nqfils(iq).gt.0) then 841 enddo 842 enddo 884 843 885 844 !write(*,*) 'vly 853: sortie' … … 888 847 END 889 848 RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq) 890 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils, ! CRisi891 & qperemin,masseqmin,ratiomin! MVals et CRisi849 USE infotrac, ONLY : nqtot,tracers, ! CRisi 850 & min_qParent,min_qMass,min_ratio ! MVals et CRisi 892 851 c 893 852 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 917 876 c --------- 918 877 c 919 INTEGER i ,ij,l,j,ii878 INTEGER ij,l 920 879 c 921 880 REAL wq(ip1jmp1,llm+1),newmasse … … 930 889 SAVE testcpu 931 890 932 REAL temps0,temps1,temps2,temps3,temps4,temps5,second 933 SAVE temps0,temps1,temps2,temps3,temps4,temps5934 REAL SSUM891 #ifdef BIDON 892 REAL temps0,temps1,second 893 SAVE temps0,temps1 935 894 936 895 DATA testcpu/.false./ 937 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./ 896 DATA temps0,temps1/0.,0./ 897 #endif 938 898 939 899 c On oriente tout dans le sens de la pression c'est a dire dans le … … 1009 969 ! CRisi: appel récursif de l'advection sur les fils. 1010 970 ! Il faut faire ça avant d'avoir mis à jour q et masse 1011 !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq) 1012 if (nqfils(iq).gt.0) then 1013 do ifils=1,nqdesc(iq) 1014 iq2=iqfils(ifils,iq) 1015 DO l=1,llm 971 !write(*,*) 'vlsplt 942: iq,nqChildren(iq)=',iq,nqChildren(iq) 972 do ifils=1,tracers(iq)%nqDescen 973 iq2=tracers(iq)%iqDescen(ifils) 974 DO l=1,llm 1016 975 DO ij=1,ip1jmp1 1017 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)1018 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)1019 !MVals: veiller a ce qu'on n'ait pas de denominateur nul1020 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)1021 if (q(ij,l,iq).gt.qperemin) then1022 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)1023 else1024 Ratio(ij,l,iq2)=ratiomin1025 endif976 !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 977 !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 978 !MVals: veiller a ce qu'on n'ait pas de denominateur nul 979 masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass) 980 if (q(ij,l,iq).gt.min_qParent) then 981 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 982 else 983 Ratio(ij,l,iq2)=min_ratio 984 endif 1026 985 enddo 1027 1028 enddo !do ifils=1,nqdesc(iq)986 enddo 987 enddo 1029 988 1030 do ifils=1,nqfils(iq) 1031 iq2=iqfils(ifils,iq) 1032 call vlz(Ratio,pente_max,masseq,wq,iq2) 1033 enddo !do ifils=1,nqfils(iq) 1034 endif !if (nqfils(iq).gt.0) then 989 do ifils=1,tracers(iq)%nqChildren 990 iq2=tracers(iq)%iqDescen(ifils) 991 call vlz(Ratio,pente_max,masseq,wq,iq2) 992 enddo 1035 993 ! end CRisi 1036 994 … … 1045 1003 1046 1004 ! retablir les fils en rapport de melange par rapport a l'air: 1047 if (nqfils(iq).gt.0) then 1048 do ifils=1,nqdesc(iq) 1049 iq2=iqfils(ifils,iq) 1050 DO l=1,llm 1005 do ifils=1,tracers(iq)%nqDescen 1006 iq2=tracers(iq)%iqDescen(ifils) 1007 DO l=1,llm 1051 1008 DO ij=1,ip1jmp1 1052 1009 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 1053 1010 enddo 1054 enddo 1055 enddo !do ifils=1,nqdesc(iq) 1056 endif !if (nqfils(iq).gt.0) then 1011 enddo 1012 enddo 1057 1013 !write(*,*) 'vlsplt 1032' 1058 1014 … … 1099 1055 real zzq(iip1,jjp1,llm) 1100 1056 1057 #ifdef isminmax 1101 1058 integer imin,jmin,lmin,ijlmin 1102 1059 integer imax,jmax,lmax,ijlmax … … 1104 1061 integer ismin,ismax 1105 1062 1106 #ifdef isminismax1107 1063 call scopy (ip1jmp1*llm,zq,1,zzq,1) 1108 1064 … … 1132 1088 #endif 1133 1089 return 1134 9999 format(a20,' q(',i3,',',i2,',',i2,')=',e12.5,e12.5)1090 c9999 format(a20,' q(',i3,',',i2,',',i2,')=',e12.5,e12.5) 1135 1091 end 1136 1092 -
LMDZ6/branches/Ocean_skin/libf/dyn3d/vlspltqs.F
r2603 r4368 4 4 SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt, 5 5 , p,pk,teta,iq ) 6 USE infotrac, ONLY: nqtot, nqdesc,iqfils6 USE infotrac, ONLY: nqtot,tracers 7 7 c 8 8 c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron … … 121 121 CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1) 122 122 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1) 123 if (nqdesc(iq).gt.0) then 124 do ifils=1,nqdesc(iq) 125 iq2=iqfils(ifils,iq) 123 do ifils=1,tracers(iq)%nqDescen 124 iq2=tracers(iq)%iqDescen(ifils) 126 125 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 127 enddo 128 endif !if (nqfils(iq).gt.0) then 126 enddo 129 127 130 128 c call minmaxq(zq,qmin,qmax,'avant vlxqs ') … … 162 160 ENDDO 163 161 ! CRisi: aussi pour les fils 164 if (nqdesc(iq).gt.0) then 165 do ifils=1,nqdesc(iq) 166 iq2=iqfils(ifils,iq) 162 do ifils=1,tracers(iq)%nqDescen 163 iq2=tracers(iq)%iqDescen(ifils) 167 164 DO l=1,llm 168 DO ij=1,ip1jmp1169 q(ij,l,iq2)=zq(ij,l,iq2)170 ENDDO171 DO ij=1,ip1jm+1,iip1165 DO ij=1,ip1jmp1 166 q(ij,l,iq2)=zq(ij,l,iq2) 167 ENDDO 168 DO ij=1,ip1jm+1,iip1 172 169 q(ij+iim,l,iq2)=q(ij,l,iq2) 173 ENDDO170 ENDDO 174 171 ENDDO 175 enddo !do ifils=1,nqdesc(iq) 176 endif ! if (nqfils(iq).gt.0) then 172 enddo 177 173 !write(*,*) 'vlspltqs 183: fin de la routine' 178 174 … … 180 176 END 181 177 SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq) 182 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils ! CRisi178 USE infotrac, ONLY : nqtot,tracers ! CRisi 183 179 184 180 c … … 483 479 ! CRisi: appel récursif de l'advection sur les fils. 484 480 ! Il faut faire ça avant d'avoir mis à jour q et masse 485 !write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq) 481 !write(*,*) 'vlspltqs 326: iq,nqChildren(iq)=',iq, 482 ! & tracers(iq)%nqChildren 486 483 487 if (nqfils(iq).gt.0) then 488 do ifils=1,nqdesc(iq) 489 iq2=iqfils(ifils,iq) 490 DO l=1,llm 484 do ifils=1,tracers(iq)%nqDescen 485 iq2=tracers(iq)%iqDescen(ifils) 486 DO l=1,llm 491 487 DO ij=iip2,ip1jm 492 493 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)494 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)488 ! On a besoin de q et masse seulement entre iip2 et ip1jm 489 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 490 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 495 491 enddo 496 enddo 497 enddo !do ifils=1,nqdesc(iq) 498 do ifils=1,nqfils(iq) 499 iq2=iqfils(ifils,iq) 500 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 501 enddo !do ifils=1,nqfils(iq) 502 endif !if (nqfils(iq).gt.0) then 492 enddo 493 enddo 494 do ifils=1,tracers(iq)%nqChildren 495 iq2=tracers(iq)%iqDescen(ifils) 496 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 497 enddo 503 498 ! end CRisi 504 499 … … 523 518 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 524 519 ! puis on boucle en longitude 525 if (nqdesc(iq).gt.0) then 526 do ifils=1,nqdesc(iq) 527 iq2=iqfils(ifils,iq) 528 DO l=1,llm 520 do ifils=1,tracers(iq)%nqDescen 521 iq2=tracers(iq)%iqDescen(ifils) 522 DO l=1,llm 529 523 DO ij=iip2+1,ip1jm 530 524 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 531 525 enddo 532 526 DO ij=iip1+iip1,ip1jm,iip1 533 q(ij-iim,l,iq2)=q(ij,l,iq2) 534 enddo ! DO ij=ijb+iip1-1,ije,iip1 535 enddo !DO l=1,llm 536 enddo !do ifils=1,nqdesc(iq) 537 endif !if (nqfils(iq).gt.0) then 527 q(ij-iim,l,iq2)=q(ij,l,iq2) 528 enddo 529 enddo 530 enddo 538 531 539 532 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) … … 544 537 END 545 538 SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq) 546 USE infotrac, ONLY : nqtot, nqfils,nqdesc,iqfils ! CRisi539 USE infotrac, ONLY : nqtot,tracers ! CRisi 547 540 c 548 541 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 794 787 ! CRisi: appel récursif de l'advection sur les fils. 795 788 ! Il faut faire ça avant d'avoir mis à jour q et masse 796 !write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq) 789 !write(*,*) 'vlyqs 689: iq,nqChildren(iq)=',iq, 790 ! & tracers(iq)%nqChildren 797 791 798 if (nqfils(iq).gt.0) then 799 do ifils=1,nqdesc(iq) 800 iq2=iqfils(ifils,iq) 801 DO l=1,llm 802 DO ij=1,ip1jmp1 803 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 804 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 792 do ifils=1,tracers(iq)%nqDescen 793 iq2=tracers(iq)%iqDescen(ifils) 794 DO l=1,llm 795 DO ij=1,ip1jmp1 796 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 797 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 805 798 enddo 806 enddo 807 enddo !do ifils=1,nqdesc(iq) 808 809 do ifils=1,nqfils(iq) 810 iq2=iqfils(ifils,iq) 811 !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2 812 call vly(Ratio,pente_max,masseq,qbyv,iq2) 813 enddo !do ifils=1,nqfils(iq) 814 endif !if (nqfils(iq).gt.0) then 799 enddo 800 enddo 801 do ifils=1,tracers(iq)%nqChildren 802 iq2=tracers(iq)%iqDescen(ifils) 803 !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2 804 call vly(Ratio,pente_max,masseq,qbyv,iq2) 805 enddo 815 806 816 807 DO l=1,llm … … 868 859 869 860 ! retablir les fils en rapport de melange par rapport a l'air: 870 if (nqdesc(iq).gt.0) then 871 do ifils=1,nqdesc(iq) 872 iq2=iqfils(ifils,iq) 873 DO l=1,llm 861 do ifils=1,tracers(iq)%nqDescen 862 iq2=tracers(iq)%iqDescen(ifils) 863 DO l=1,llm 874 864 DO ij=1,ip1jmp1 875 865 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 876 866 enddo 877 enddo 878 enddo !do ifils=1,nqdesc(iq) 879 endif !if (nqfils(iq).gt.0) then 867 enddo 868 enddo 880 869 !write(*,*) 'vly 879' 881 870
Note: See TracChangeset
for help on using the changeset viewer.
