- Timestamp:
- Jan 15, 2016, 8:27:16 AM (9 years ago)
- Location:
- trunk
- Files:
-
- 4 added
- 1 deleted
- 13 edited
- 6 moved
-
DOC/chantiers/commit_importants.log (modified) (1 diff)
-
LMDZ.COMMON/libf/dyn3d/advtrac.F90 (modified) (5 diffs)
-
LMDZ.COMMON/libf/dyn3d/caladvtrac.F (modified) (4 diffs)
-
LMDZ.COMMON/libf/dyn3d/check_isotopes.F (added)
-
LMDZ.COMMON/libf/dyn3d/guide_mod.F90 (modified) (2 diffs)
-
LMDZ.COMMON/libf/dyn3d/leapfrog.F (modified) (12 diffs)
-
LMDZ.COMMON/libf/dyn3d/qminimum.F (modified) (7 diffs)
-
LMDZ.COMMON/libf/dyn3d/vlsplt.F (modified) (37 diffs)
-
LMDZ.COMMON/libf/dyn3d/vlspltqs.F (modified) (28 diffs)
-
LMDZ.COMMON/libf/dyn3d_common/covcont.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/covcont.F) (1 diff)
-
LMDZ.COMMON/libf/dyn3d_common/dynetat0.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/dynetat0.F) (8 diffs)
-
LMDZ.COMMON/libf/dyn3d_common/dynredem.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/dynredem.F) (3 diffs)
-
LMDZ.COMMON/libf/dyn3d_common/dynredem_mod.F90 (added)
-
LMDZ.COMMON/libf/dyn3d_common/grid_atob.F (deleted)
-
LMDZ.COMMON/libf/dyn3d_common/infotrac.F90 (modified) (18 diffs)
-
LMDZ.COMMON/libf/dyn3d_common/iniacademic.F90 (modified) (2 diffs)
-
LMDZ.COMMON/libf/dyn3d_common/iso_verif_dyn.F (added)
-
LMDZ.COMMON/libf/dyn3d_common/misc_mod.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/misc/misc_mod.F90)
-
LMDZ.COMMON/libf/dyn3d_common/pression.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3d_common/pression.F) (1 diff)
-
LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F) (3 diffs)
-
LMDZ.COMMON/libf/dyn3dpar/integrd_p.F (modified) (1 diff)
-
LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F (modified) (2 diffs)
-
LMDZ.COMMON/libf/dynlonlat_phylonlat/grid_atob_m.F90 (added)
-
LMDZ.COMMON/libf/misc/wxios.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified trunk/DOC/chantiers/commit_importants.log ¶
r1441 r1508 1558 1558 - ran1.F, sort.F, minmax.F, minmax2.F, juldate.F moved over from dyn3d_common 1559 1559 1560 ********************** 1561 **** commit_v1508 **** 1562 ********************** 1563 Ehouarn: Updates in the dynamics (seq and //) to keep up with updates 1564 in LMDZ5 (up to LMDZ5 trunk, rev 2325): 1565 IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar 1566 as in LMDZ5 these modifications were done in dyn3dmem. 1567 Related LMDZ5 revisions are r2270 and r2281 1568 * in dynlonlat_phylonlat: 1569 - add module "grid_atob_m.F90" (a regridding utility 1570 so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F) 1571 1572 * in misc: 1573 - follow up updates on wxios.F (add missing_val module variable) 1574 1575 * in dyn3d_common: 1576 - pression.F => pression.F90 1577 - misc_mod.F90: moved from misc to dyn3d_common 1578 - added new iso_verif_dyn.F 1579 - covcont.F => covcont.F90 1580 - infotrac.F90 : add handling of isotopes (reading of corresponding 1581 traceur.def for planets not implemented) 1582 - dynetat0.F => dynetat0.F90 with some code factorization 1583 - dynredem.F => dynredem.F90 with some code factorization 1584 - added dynredem_mod.F90: routines used by dynredem 1585 - iniacademic.F90 : added isotopes-related initialization for Earth case 1586 1587 * in dyn3d: 1588 - added check_isotopes.F 1589 - modified (isotopes) advtrac.F90, caladvtrac.F 1590 - guide_mod.F90: ported updates 1591 - leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers) 1592 - qminimium.F : adaptations for isotopes (copied over, except that #include 1593 comvert.h is not needed). 1594 - vlsplt.F: adaptations for isotopes (copied over, except than #include 1595 logic.h, comvert.h not needed, and replace "include comconst.h" with 1596 use comconst_mod, ONLY: pi) 1597 - vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T) 1598 1599 * in dyn3dpar: 1600 - leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call 1601 integrd_p with nqtot tracers (only important for Earth) 1602 - dynredem_p.F => dynredem_p.F90 and some code factorization 1603 - and no isotopes-relates changes in dyn3dpar (since these changes 1604 have been made in LMDZ5 dyn3dmem). -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90 ¶
r1441 r1508 9 9 ! M.A Filiberti (04/2002) 10 10 ! 11 USE infotrac, ONLY: nqtot, iadv 11 USE infotrac, ONLY: nqtot, iadv, nqperes, ok_iso_verif 12 12 USE control_mod, ONLY: iapp_tracvl, day_step 13 13 USE comconst_mod, ONLY: dtvr … … 218 218 ! Appel des sous programmes d'advection 219 219 !----------------------------------------------------------- 220 do iq=1,nqtot 220 if (ok_iso_verif) then 221 write(*,*) 'advtrac 227' 222 call check_isotopes_seq(q,ip1jmp1,'advtrac 162') 223 endif !if (ok_iso_verif) then 224 225 do iq=1,nqperes 221 226 ! call clock(t_initial) 222 227 if(iadv(iq) == 0) cycle … … 225 230 ! ---------------------------------------------------------------- 226 231 if(iadv(iq).eq.10) THEN 227 call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr) 232 ! CRisi: on fait passer tout q pour avoir acces aux fils 233 234 !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:) 235 call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) 228 236 229 237 ! ---------------------------------------------------------------- … … 233 241 else if(iadv(iq).eq.14) then 234 242 ! 235 CALL vlspltqs( q(1,1,1), 2., massem, wg , & 236 pbarug,pbarvg,dtvr,p,pk,teta ) 243 !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:) 244 CALL vlspltqs( q, 2., massem, wg , & 245 pbarug,pbarvg,dtvr,p,pk,teta,iq) 246 237 247 ! ---------------------------------------------------------------- 238 248 ! Schema de Frederic Hourdin … … 383 393 end DO 384 394 395 if (ok_iso_verif) then 396 write(*,*) 'advtrac 402' 397 call check_isotopes_seq(q,ip1jmp1,'advtrac 397') 398 endif !if (ok_iso_verif) then 385 399 386 400 !------------------------------------------------------------------ -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d/caladvtrac.F ¶
r1422 r1508 53 53 if (planet_type.eq."earth") then 54 54 C initialisation 55 dq(:,:,1:2)=q(:,:,1:2) 55 ! CRisi: il faut gérer tous les traceurs si on veut pouvoir faire des 56 ! isotopes 57 ! dq(:,:,1:2)=q(:,:,1:2) 58 dq(:,:,1:nqtot)=q(:,:,1:nqtot) 56 59 57 60 c test des valeurs minmax … … 82 85 ENDDO 83 86 84 CALL qminimum( q, 2, finmasse ) 87 !write(*,*) 'caladvtrac 87' 88 CALL qminimum( q, nqtot, finmasse ) 89 !write(*,*) 'caladvtrac 89' 85 90 86 91 CALL SCOPY ( ip1jmp1*llm, masse, 1, finmasse, 1 ) … … 92 97 dtvrtrac = iapp_tracvl * dtvr 93 98 c 94 DO iq = 1 , 299 DO iq = 1 , nqtot 95 100 DO l = 1 , llm 96 101 DO ij = 1,ip1jmp1 … … 105 110 if (planet_type.eq."earth") then 106 111 ! Earth-specific treatment for the first 2 tracers (water) 107 dq(:,:,1: 2)=0.112 dq(:,:,1:nqtot)=0. 108 113 endif ! of if (planet_type.eq."earth") 109 114 ENDIF ! of IF( iapptrac.EQ.iapp_tracvl ) -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90 ¶
r1441 r1508 79 79 ! Lecture des parametres: 80 80 ! --------------------------------------------- 81 call ini_getparam("nudging_parameters_out.txt") 81 82 ! Variables guidees 82 83 CALL getpar('guide_u',.true.,guide_u,'guidage de u') … … 144 145 CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S') 145 146 CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P') 147 148 call fin_getparam 146 149 147 150 ! --------------------------------------------- -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F ¶
r1422 r1508 12 12 use IOIPSL 13 13 #endif 14 USE infotrac, ONLY: nqtot 14 USE infotrac, ONLY: nqtot,ok_iso_verif 15 15 USE guide_mod, ONLY : guide_main 16 16 USE write_field, ONLY: writefield … … 305 305 jH_cur = jH_cur - int(jH_cur) 306 306 307 if (ok_iso_verif) then 308 call check_isotopes_seq(q,ip1jmp1,'leapfrog 321') 309 endif !if (ok_iso_verif) then 307 310 308 311 #ifdef CPP_IOIPSL … … 337 340 ! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 ) 338 341 ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 ) 342 343 if (ok_iso_verif) then 344 call check_isotopes_seq(q,ip1jmp1,'leapfrog 400') 345 endif !if (ok_iso_verif) then 339 346 340 347 2 CONTINUE ! Matsuno backward or leapfrog step begins here … … 381 388 #endif 382 389 390 if (ok_iso_verif) then 391 call check_isotopes_seq(q,ip1jmp1,'leapfrog 589') 392 endif !if (ok_iso_verif) then 393 383 394 c----------------------------------------------------------------------- 384 395 c calcul des tendances dynamiques: … … 419 430 c ------------------------------------------------------------- 420 431 432 if (ok_iso_verif) then 433 call check_isotopes_seq(q,ip1jmp1, 434 & 'leapfrog 686: avant caladvtrac') 435 endif !if (ok_iso_verif) then 436 421 437 IF( forward. OR . leapf ) THEN 422 438 ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step … … 443 459 c ---------------------------------- 444 460 445 446 CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 , 461 if (ok_iso_verif) then 462 write(*,*) 'leapfrog 720' 463 call check_isotopes_seq(q,ip1jmp1,'leapfrog 756') 464 endif !if (ok_iso_verif) then 465 466 CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 , 447 467 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ) 448 468 ! $ finvmaold ) 469 470 if (ok_iso_verif) then 471 write(*,*) 'leapfrog 724' 472 call check_isotopes_seq(q,ip1jmp1,'leapfrog 762') 473 endif !if (ok_iso_verif) then 449 474 450 475 IF ((planet_type.eq."titan").and.(tidal)) then … … 519 544 IF (ip_ebil_dyn.ge.1 ) THEN 520 545 ztit='bil dyn' 521 ! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!546 ! Ehouarn: be careful, diagedyn is Earth-specific! 522 547 IF (planet_type.eq."earth") THEN 523 548 CALL diagedyn(ztit,2,1,1,dtphys … … 619 644 CALL massdair(p,masse) 620 645 646 if (ok_iso_verif) then 647 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196') 648 endif !if (ok_iso_verif) then 649 621 650 c----------------------------------------------------------------------- 622 651 c dissipation horizontale et verticale des petites echelles: … … 717 746 c preparation du pas d'integration suivant ...... 718 747 748 if (ok_iso_verif) then 749 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509') 750 endif !if (ok_iso_verif) then 751 719 752 IF ( .NOT.purmats ) THEN 720 753 c ........................................................ … … 776 809 777 810 ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin)) 811 812 if (ok_iso_verif) then 813 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584') 814 endif !if (ok_iso_verif) then 778 815 779 816 c----------------------------------------------------------------------- … … 858 895 ELSE ! of IF (.not.purmats) 859 896 897 if (ok_iso_verif) then 898 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664') 899 endif !if (ok_iso_verif) then 900 860 901 c ........................................................ 861 902 c .............. schema matsuno ............... … … 880 921 881 922 ELSE ! of IF(forward) i.e. backward step 923 924 if (ok_iso_verif) then 925 call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698') 926 endif !if (ok_iso_verif) then 882 927 883 928 IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d/qminimum.F ¶
r1422 r1508 2 2 ! $Header$ 3 3 ! 4 SUBROUTINE qminimum( q,nq ,deltap )4 SUBROUTINE qminimum( q,nqtot,deltap ) 5 5 6 USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif 6 7 IMPLICIT none 7 8 c … … 12 13 #include "paramet.h" 13 14 c 14 INTEGER nq 15 REAL q(ip1jmp1,llm,nq ), deltap(ip1jmp1,llm)15 INTEGER nqtot 16 REAL q(ip1jmp1,llm,nqtot), deltap(ip1jmp1,llm) 16 17 c 17 18 INTEGER iq_vap, iq_liq … … 29 30 INTEGER i, k, iq 30 31 REAL zx_defau, zx_abc, zx_pump(ip1jmp1), pompe 32 33 real zx_defau_diag(ip1jmp1,llm,2) 34 real q_follow(ip1jmp1,llm,2) 31 35 c 32 36 REAL SSUM … … 35 39 SAVE imprim 36 40 DATA imprim /0/ 41 !INTEGER ijb,ije 42 !INTEGER Index_pump(ij_end-ij_begin+1) 43 !INTEGER nb_pump 44 INTEGER ixt 37 45 c 38 46 c Quand l'eau liquide est trop petite (ou negative), on prend … … 40 48 c (sans changer la temperature !) 41 49 c 50 51 if (ok_iso_verif) then 52 call check_isotopes_seq(q,ip1jmp1,'qminimum 52') 53 endif !if (ok_iso_verif) then 54 55 zx_defau_diag(:,:,:)=0.0 56 q_follow(:,:,1:2)=q(:,:,1:2) 42 57 DO 1000 k = 1, llm 43 58 DO 1040 i = 1, ip1jmp1 44 59 if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then 60 61 if (ok_isotopes) then 62 zx_defau_diag(i,k,iq_liq)=AMAX1 63 : ( seuil_liq - q(i,k,iq_liq), 0.0 ) 64 endif !if (ok_isotopes) then 65 45 66 q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq 46 67 q(i,k,iq_liq) = seuil_liq … … 58 79 DO i = 1, ip1jmp1 59 80 if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then 81 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 85 60 86 q(i,k-1,iq) = q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) * 61 87 & deltap(i,k) / deltap(i,k-1) … … 82 108 ENDDO 83 109 ENDIF 110 111 !write(*,*) 'qminimum 128' 112 if (ok_isotopes) then 113 ! CRisi: traiter de même les traceurs d'eau 114 ! Mais il faut les prendre à l'envers pour essayer de conserver la 115 ! masse. 116 ! 1) pompage dans le sol 117 ! On suppose que ce pompage se fait sans isotopes -> on ne modifie 118 ! rien ici et on croise les doigts pour que ça ne soit pas trop 119 ! génant 120 DO i = 1,ip1jmp1 121 if (zx_pump(i).gt.0.0) then 122 q_follow(i,1,iq_vap)=q_follow(i,1,iq_vap)+zx_pump(i) 123 endif !if (zx_pump(i).gt.0.0) then 124 enddo !DO i = 1,ip1jmp1 125 126 ! 2) transfert de vap vers les couches plus hautes 127 !write(*,*) 'qminimum 139' 128 do k=2,llm 129 DO i = 1,ip1jmp1 130 if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 131 ! on ajoute la vapeur en k 132 do ixt=1,ntraciso 133 q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap)) 134 : +zx_defau_diag(i,k,iq_vap) 135 : *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap) 136 137 ! 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 : -zx_defau_diag(i,k,iq_vap) 140 : *deltap(i,k)/deltap(i,k-1) 141 : *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap) 142 143 enddo !do ixt=1,niso 144 q_follow(i,k,iq_vap)= q_follow(i,k,iq_vap) 145 : +zx_defau_diag(i,k,iq_vap) 146 q_follow(i,k-1,iq_vap)= q_follow(i,k-1,iq_vap) 147 : -zx_defau_diag(i,k,iq_vap) 148 : *deltap(i,k)/deltap(i,k-1) 149 endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 150 enddo !DO i = 1, ip1jmp1 151 enddo !do k=2,llm 152 153 if (ok_iso_verif) then 154 call check_isotopes_seq(q,ip1jmp1,'qminimum 168') 155 endif !if (ok_iso_verif) then 156 157 158 ! 3) transfert d'eau de la vapeur au liquide 159 !write(*,*) 'qminimum 164' 160 do k=1,llm 161 DO i = 1,ip1jmp1 162 if (zx_defau_diag(i,k,iq_liq).gt.0.0) then 163 164 ! on ajoute eau liquide en k en k 165 do ixt=1,ntraciso 166 q(i,k,iqiso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq)) 167 : +zx_defau_diag(i,k,iq_liq) 168 : *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap) 169 ! et on la retranche à la vapeur en k 170 q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap)) 171 : -zx_defau_diag(i,k,iq_liq) 172 : *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap) 173 enddo !do ixt=1,niso 174 q_follow(i,k,iq_liq)= q_follow(i,k,iq_liq) 175 : +zx_defau_diag(i,k,iq_liq) 176 q_follow(i,k,iq_vap)= q_follow(i,k,iq_vap) 177 : -zx_defau_diag(i,k,iq_liq) 178 endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then 179 enddo !DO i = 1, ip1jmp1 180 enddo !do k=2,llm 181 182 if (ok_iso_verif) then 183 call check_isotopes_seq(q,ip1jmp1,'qminimum 197') 184 endif !if (ok_iso_verif) then 185 186 endif !if (ok_isotopes) then 187 !write(*,*) 'qminimum 188' 188 84 189 c 85 190 RETURN -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d/vlsplt.F ¶
r1422 r1508 1 1 c 2 c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z lguez $ 3 c 4 5 SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt) 2 c $Id: vlsplt.F 2286 2015-05-20 13:27:07Z emillour $ 3 c 4 5 SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq) 6 USE infotrac, ONLY: nqtot,nqdesc,iqfils 6 7 c 7 8 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 29 30 c REAL masse(iip1,jjp1,llm),pente_max 30 31 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm) 31 REAL q(ip1jmp1,llm )32 REAL q(ip1jmp1,llm,nqtot) 32 33 c REAL q(iip1,jjp1,llm) 33 34 REAL w(ip1jmp1,llm),pdt 35 INTEGER iq ! CRisi 34 36 c 35 37 c Local … … 39 41 INTEGER ijlqmin,iqmin,jqmin,lqmin 40 42 c 41 REAL zm(ip1jmp1,llm ),newmasse43 REAL zm(ip1jmp1,llm,nqtot),newmasse 42 44 REAL mu(ip1jmp1,llm) 43 45 REAL mv(ip1jm,llm) 44 46 REAL mw(ip1jmp1,llm+1) 45 REAL zq(ip1jmp1,llm ),zz47 REAL zq(ip1jmp1,llm,nqtot),zz 46 48 REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm) 47 49 REAL second,temps0,temps1,temps2,temps3 … … 52 54 SAVE temps1,temps2,temps3 53 55 INTEGER iminn,imaxx 56 INTEGER ifils,iq2 ! CRisi 54 57 55 58 REAL qmin,qmax … … 76 79 mw(ij,llm+1)=0. 77 80 ENDDO 78 79 CALL SCOPY(ijp1llm,q,1,zq,1) 80 CALL SCOPY(ijp1llm,masse,1,zm,1) 81 82 CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1) 83 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1) 84 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 81 91 82 92 cprint*,'Entree vlx1' 83 93 c call minmaxq(zq,qmin,qmax,'avant vlx ') 84 call vlx(zq,pente_max,zm,mu )94 call vlx(zq,pente_max,zm,mu,iq) 85 95 cprint*,'Sortie vlx1' 86 96 c call minmaxq(zq,qmin,qmax,'apres vlx1 ') 87 97 88 98 c print*,'Entree vly1' 89 call vly(zq,pente_max,zm,mv) 99 100 call vly(zq,pente_max,zm,mv,iq) 90 101 c call minmaxq(zq,qmin,qmax,'apres vly1 ') 91 102 cprint*,'Sortie vly1' 92 call vlz(zq,pente_max,zm,mw )103 call vlz(zq,pente_max,zm,mw,iq) 93 104 c call minmaxq(zq,qmin,qmax,'apres vlz ') 94 105 95 106 96 call vly(zq,pente_max,zm,mv )107 call vly(zq,pente_max,zm,mv,iq) 97 108 c call minmaxq(zq,qmin,qmax,'apres vly ') 98 109 99 110 100 call vlx(zq,pente_max,zm,mu )111 call vlx(zq,pente_max,zm,mu,iq) 101 112 c call minmaxq(zq,qmin,qmax,'apres vlx2 ') 102 113 … … 104 115 DO l=1,llm 105 116 DO ij=1,ip1jmp1 106 q(ij,l )=zq(ij,l)117 q(ij,l,iq)=zq(ij,l,iq) 107 118 ENDDO 108 119 DO ij=1,ip1jm+1,iip1 109 q(ij+iim,l)=q(ij,l) 110 ENDDO 111 ENDDO 120 q(ij+iim,l,iq)=q(ij,l,iq) 121 ENDDO 122 ENDDO 123 ! CRisi: aussi pour les fils 124 if (nqdesc(iq).gt.0) then 125 do ifils=1,nqdesc(iq) 126 iq2=iqfils(ifils,iq) 127 DO l=1,llm 128 DO ij=1,ip1jmp1 129 q(ij,l,iq2)=zq(ij,l,iq2) 130 ENDDO 131 DO ij=1,ip1jm+1,iip1 132 q(ij+iim,l,iq2)=q(ij,l,iq2) 133 ENDDO 134 ENDDO 135 enddo !do ifils=1,nqdesc(iq) 136 endif ! if (nqdesc(iq).gt.0) then 112 137 113 138 RETURN 114 139 END 115 SUBROUTINE vlx(q,pente_max,masse,u_m) 140 RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq) 141 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 116 142 117 143 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 133 159 c Arguments: 134 160 c ---------- 135 REAL masse(ip1jmp1,llm ),pente_max161 REAL masse(ip1jmp1,llm,nqtot),pente_max 136 162 REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm) 137 REAL q(ip1jmp1,llm )163 REAL q(ip1jmp1,llm,nqtot) 138 164 REAL w(ip1jmp1,llm) 165 INTEGER iq ! CRisi 139 166 c 140 167 c Local … … 149 176 REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm) 150 177 REAL u_mq(ip1jmp1,llm) 178 179 ! CRisi 180 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) 181 INTEGER ifils,iq2 ! CRisi 151 182 152 183 Logical extremum,first,testcpu … … 182 213 DO l = 1, llm 183 214 DO ij=iip2,ip1jm-1 184 dxqu(ij)=q(ij+1,l )-q(ij,l)215 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 185 216 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0' 186 c sigu(ij)=u_m(ij,l)/masse(ij,l )217 c sigu(ij)=u_m(ij,l)/masse(ij,l,iq) 187 218 ENDDO 188 219 DO ij=iip1+iip1,ip1jm,iip1 … … 237 268 DO l = 1, llm 238 269 DO ij=iip2,ip1jm-1 239 dxqu(ij)=q(ij+1,l )-q(ij,l)270 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 240 271 ENDDO 241 272 DO ij=iip1+iip1,ip1jm,iip1 … … 279 310 DO l=1,llm 280 311 DO ij=iip2,ip1jm-1 281 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l ),282 , 1.+u_m(ij,l)/masse(ij+1,l ),312 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq), 313 , 1.+u_m(ij,l)/masse(ij+1,l,iq), 283 314 , u_m(ij,l)) 284 315 zdum(ij,l)=0.5*zdum(ij,l) 285 316 u_mq(ij,l)=cvmgp( 286 , q(ij,l )+zdum(ij,l)*dxq(ij,l),287 , q(ij+1,l )-zdum(ij,l)*dxq(ij+1,l),317 , q(ij,l,iq)+zdum(ij,l)*dxq(ij,l), 318 , q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l), 288 319 , u_m(ij,l)) 289 320 u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l) … … 297 328 DO l=1,llm 298 329 DO ij=iip2,ip1jm-1 299 c print*,'masse(',ij,')=',masse(ij,l )330 c print*,'masse(',ij,')=',masse(ij,l,iq) 300 331 IF (u_m(ij,l).gt.0.) THEN 301 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l )302 u_mq(ij,l)=u_m(ij,l)*(q(ij,l )+0.5*zdum(ij,l)*dxq(ij,l))332 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq) 333 u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l)) 303 334 ELSE 304 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l) 305 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l)) 335 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq) 336 u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq) 337 & -0.5*zdum(ij,l)*dxq(ij+1,l)) 306 338 ENDIF 307 339 ENDDO … … 373 405 i=ijq-(j-1)*iip1 374 406 c accumulation pour les mailles completements advectees 375 do while(zu_m.gt.masse(ijq,l)) 376 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l) 377 zu_m=zu_m-masse(ijq,l) 407 do while(zu_m.gt.masse(ijq,l,iq)) 408 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq) 409 & *masse(ijq,l,iq) 410 zu_m=zu_m-masse(ijq,l,iq) 378 411 i=mod(i-2+iim,iim)+1 379 412 ijq=(j-1)*iip1+i … … 381 414 c ajout de la maille non completement advectee 382 415 u_mq(ij,l)=u_mq(ij,l)+zu_m* 383 & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l)) 416 & (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq)) 417 & *dxq(ijq,l)) 384 418 ELSE 385 419 ijq=ij+1 386 420 i=ijq-(j-1)*iip1 387 421 c accumulation pour les mailles completements advectees 388 do while(-zu_m.gt.masse(ijq,l)) 389 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l) 390 zu_m=zu_m+masse(ijq,l) 422 do while(-zu_m.gt.masse(ijq,l,iq)) 423 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) 424 & *masse(ijq,l,iq) 425 zu_m=zu_m+masse(ijq,l,iq) 391 426 i=mod(i,iim)+1 392 427 ijq=(j-1)*iip1+i 393 428 ENDDO 394 429 c ajout de la maille non completement advectee 395 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l )-396 & 0.5*(1.+zu_m/masse(ijq,l ))*dxq(ijq,l))430 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- 431 & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 397 432 ENDIF 398 433 ENDDO … … 411 446 ENDDO 412 447 448 ! CRisi: appel récursif de l'advection sur les fils. 449 ! Il faut faire ça avant d'avoir mis à jour q et masse 450 !write(*,*) 'vlsplt 326: iq,nqfils(iq)=',iq,nqfils(iq) 451 452 if (nqdesc(iq).gt.0) then 453 do ifils=1,nqdesc(iq) 454 iq2=iqfils(ifils,iq) 455 DO l=1,llm 456 DO ij=iip2,ip1jm 457 ! On a besoin de q et masse seulement entre iip2 et ip1jm 458 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 459 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 460 enddo 461 enddo 462 enddo !do ifils=1,nqdesc(iq) 463 do ifils=1,nqfils(iq) 464 iq2=iqfils(ifils,iq) 465 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 466 enddo !do ifils=1,nqfils(iq) 467 endif !if (nqfils(iq).gt.0) then 468 ! end CRisi 469 413 470 414 471 c calcul des tENDances … … 416 473 DO l=1,llm 417 474 DO ij=iip2+1,ip1jm 418 new_m=masse(ij,l )+u_m(ij-1,l)-u_m(ij,l)419 q(ij,l )=(q(ij,l)*masse(ij,l)+475 new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l) 476 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ 420 477 & u_mq(ij-1,l)-u_mq(ij,l)) 421 478 & /new_m 422 masse(ij,l )=new_m479 masse(ij,l,iq)=new_m 423 480 ENDDO 424 481 c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous) 425 482 DO ij=iip1+iip1,ip1jm,iip1 426 q(ij-iim,l)=q(ij,l) 427 masse(ij-iim,l)=masse(ij,l) 428 ENDDO 429 ENDDO 483 q(ij-iim,l,iq)=q(ij,l,iq) 484 masse(ij-iim,l,iq)=masse(ij,l,iq) 485 ENDDO 486 ENDDO 487 488 ! retablir les fils en rapport de melange par rapport a l'air: 489 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 490 ! puis on boucle en longitude 491 if (nqdesc(iq).gt.0) then 492 do ifils=1,nqdesc(iq) 493 iq2=iqfils(ifils,iq) 494 DO l=1,llm 495 DO ij=iip2+1,ip1jm 496 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 497 enddo 498 DO ij=iip1+iip1,ip1jm,iip1 499 q(ij-iim,l,iq2)=q(ij,l,iq2) 500 enddo ! DO ij=ijb+iip1-1,ije,iip1 501 enddo !DO l=1,llm 502 enddo !do ifils=1,nqdesc(iq) 503 endif !if (nqfils(iq).gt.0) then 504 430 505 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) 431 506 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1) … … 434 509 RETURN 435 510 END 436 SUBROUTINE vly(q,pente_max,masse,masse_adv_v) 511 RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq) 512 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 513 USE comconst_mod, ONLY: pi 437 514 c 438 515 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 446 523 c 447 524 c -------------------------------------------------------------------- 448 USE comconst_mod, ONLY: pi449 450 525 IMPLICIT NONE 451 526 c … … 457 532 c Arguments: 458 533 c ---------- 459 REAL masse(ip1jmp1,llm ),pente_max534 REAL masse(ip1jmp1,llm,nqtot),pente_max 460 535 REAL masse_adv_v( ip1jm,llm) 461 REAL q(ip1jmp1,llm), dq( ip1jmp1,llm) 536 REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm) 537 INTEGER iq ! CRisi 462 538 c 463 539 c Local … … 484 560 SAVE sinlon,coslon,sinlondlon,coslondlon 485 561 SAVE airej2,airejjm 562 563 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi 564 INTEGER ifils,iq2 ! CRisi 565 486 566 c 487 567 c … … 490 570 DATA first,testcpu/.true.,.false./ 491 571 DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./ 572 573 !write(*,*) 'vly 578: entree, iq=',iq 492 574 493 575 IF(first) THEN … … 522 604 523 605 DO i = 1, iim 524 airescb(i) = aire(i+ iip1) * q(i+ iip1,l )525 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l )606 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq) 607 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq) 526 608 ENDDO 527 609 qpns = SSUM( iim, airescb ,1 ) / airej2 … … 531 613 532 614 DO ij=1,ip1jm 533 dyqv(ij)=q(ij,l )-q(ij+iip1,l)615 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq) 534 616 adyqv(ij)=abs(dyqv(ij)) 535 617 ENDDO … … 546 628 547 629 DO ij=1,iip1 548 dyq(ij,l)=qpns-q(ij+iip1,l )549 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l )-qpsn630 dyq(ij,l)=qpns-q(ij+iip1,l,iq) 631 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn 550 632 ENDDO 551 633 … … 668 750 ENDDO 669 751 752 !write(*,*) 'vly 756' 670 753 DO l=1,llm 671 754 DO ij=1,ip1jm 672 755 IF(masse_adv_v(ij,l).gt.0) THEN 673 qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)* 674 , 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)) 756 qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)* 757 , 0.5*(1.-masse_adv_v(ij,l) 758 , /masse(ij+iip1,l,iq)) 675 759 ELSE 676 qbyv(ij,l)=q(ij,l)-dyq(ij,l)* 677 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) 760 qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)* 761 , 0.5*(1.+masse_adv_v(ij,l) 762 , /masse(ij,l,iq)) 678 763 ENDIF 679 764 qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l) … … 681 766 ENDDO 682 767 768 ! CRisi: appel récursif de l'advection sur les fils. 769 ! Il faut faire ça avant d'avoir mis à jour q et masse 770 !write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq) 771 772 if (nqfils(iq).gt.0) then 773 do ifils=1,nqdesc(iq) 774 iq2=iqfils(ifils,iq) 775 DO l=1,llm 776 DO ij=1,ip1jmp1 777 ! attention, chaque fils doit avoir son masseq, sinon, le 1er 778 ! fils ecrase le masseq de ses freres. 779 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 780 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 781 enddo 782 enddo 783 enddo !do ifils=1,nqdesc(iq) 784 785 do ifils=1,nqfils(iq) 786 iq2=iqfils(ifils,iq) 787 call vly(Ratio,pente_max,masseq,qbyv,iq2) 788 enddo !do ifils=1,nqfils(iq) 789 endif !if (nqfils(iq).gt.0) then 683 790 684 791 DO l=1,llm 685 792 DO ij=iip2,ip1jm 686 newmasse=masse(ij,l )793 newmasse=masse(ij,l,iq) 687 794 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 688 q(ij,l )=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))689 & /newmasse690 masse(ij,l )=newmasse795 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l) 796 & -qbyv(ij-iip1,l))/newmasse 797 masse(ij,l,iq)=newmasse 691 798 ENDDO 692 799 c.-. ancienne version … … 696 803 convpn=SSUM(iim,qbyv(1,l),1) 697 804 convmpn=ssum(iim,masse_adv_v(1,l),1) 698 massepn=ssum(iim,masse(1,l ),1)805 massepn=ssum(iim,masse(1,l,iq),1) 699 806 qpn=0. 700 807 do ij=1,iim 701 qpn=qpn+masse(ij,l )*q(ij,l)808 qpn=qpn+masse(ij,l,iq)*q(ij,l,iq) 702 809 enddo 703 810 qpn=(qpn+convpn)/(massepn+convmpn) 704 811 do ij=1,iip1 705 q(ij,l )=qpn812 q(ij,l,iq)=qpn 706 813 enddo 707 814 … … 711 818 convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) 712 819 convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1) 713 masseps=ssum(iim, masse(ip1jm+1,l ),1)820 masseps=ssum(iim, masse(ip1jm+1,l,iq),1) 714 821 qps=0. 715 822 do ij = ip1jm+1,ip1jmp1-1 716 qps=qps+masse(ij,l )*q(ij,l)823 qps=qps+masse(ij,l,iq)*q(ij,l,iq) 717 824 enddo 718 825 qps=(qps+convps)/(masseps+convmps) 719 826 do ij=ip1jm+1,ip1jmp1 720 q(ij,l )=qps827 q(ij,l,iq)=qps 721 828 enddo 722 829 … … 732 839 c DO ij = 1,iip1 733 840 c q(ij,l)=newq 734 c masse(ij,l )=newmasse*aire(ij)841 c masse(ij,l,iq)=newmasse*aire(ij) 735 842 c ENDDO 736 843 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) … … 742 849 c DO ij = ip1jm+1,ip1jmp1 743 850 c q(ij,l)=newq 744 c masse(ij,l )=newmasse*aire(ij)851 c masse(ij,l,iq)=newmasse*aire(ij) 745 852 c ENDDO 746 853 c._. fin nouvelle version 747 854 ENDDO 855 856 ! retablir les fils en rapport de melange par rapport a l'air: 857 if (nqfils(iq).gt.0) then 858 do ifils=1,nqdesc(iq) 859 iq2=iqfils(ifils,iq) 860 DO l=1,llm 861 DO ij=1,ip1jmp1 862 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 863 enddo 864 enddo 865 enddo !do ifils=1,nqdesc(iq) 866 endif !if (nqfils(iq).gt.0) then 867 868 !write(*,*) 'vly 853: sortie' 748 869 749 870 RETURN 750 871 END 751 SUBROUTINE vlz(q,pente_max,masse,w) 872 RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq) 873 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 752 874 c 753 875 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 769 891 c Arguments: 770 892 c ---------- 771 REAL masse(ip1jmp1,llm ),pente_max772 REAL q(ip1jmp1,llm )893 REAL masse(ip1jmp1,llm,nqtot),pente_max 894 REAL q(ip1jmp1,llm,nqtot) 773 895 REAL w(ip1jmp1,llm+1) 896 INTEGER iq 774 897 c 775 898 c Local … … 782 905 REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax 783 906 REAL sigw 907 908 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi 909 INTEGER ifils,iq2 ! CRisi 784 910 785 911 LOGICAL testcpu … … 795 921 c On oriente tout dans le sens de la pression c'est a dire dans le 796 922 c sens de W 923 924 !write(*,*) 'vlz 923: entree' 797 925 798 926 #ifdef BIDON … … 803 931 DO l=2,llm 804 932 DO ij=1,ip1jmp1 805 dzqw(ij,l)=q(ij,l-1 )-q(ij,l)933 dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq) 806 934 adzqw(ij,l)=abs(dzqw(ij,l)) 807 935 ENDDO … … 825 953 ENDDO 826 954 955 !write(*,*) 'vlz 954' 827 956 DO ij=1,ip1jmp1 828 957 dzq(ij,1)=0. … … 841 970 c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq 842 971 972 !write(*,*) 'vlz 969' 843 973 DO l = 1,llm-1 844 974 do ij = 1,ip1jmp1 845 975 IF(w(ij,l+1).gt.0.) THEN 846 sigw=w(ij,l+1)/masse(ij,l+1) 847 wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 976 sigw=w(ij,l+1)/masse(ij,l+1,iq) 977 wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq) 978 & +0.5*(1.-sigw)*dzq(ij,l+1)) 848 979 ELSE 849 sigw=w(ij,l+1)/masse(ij,l )850 wq(ij,l+1)=w(ij,l+1)*(q(ij,l )-0.5*(1.+sigw)*dzq(ij,l))980 sigw=w(ij,l+1)/masse(ij,l,iq) 981 wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l)) 851 982 ENDIF 852 983 ENDDO … … 858 989 ENDDO 859 990 991 ! CRisi: appel récursif de l'advection sur les fils. 992 ! Il faut faire ça avant d'avoir mis à jour q et masse 993 !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq) 994 if (nqfils(iq).gt.0) then 995 do ifils=1,nqdesc(iq) 996 iq2=iqfils(ifils,iq) 997 DO l=1,llm 998 DO ij=1,ip1jmp1 999 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 1000 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 1001 enddo 1002 enddo 1003 enddo !do ifils=1,nqdesc(iq) 1004 1005 do ifils=1,nqfils(iq) 1006 iq2=iqfils(ifils,iq) 1007 call vlz(Ratio,pente_max,masseq,wq,iq2) 1008 enddo !do ifils=1,nqfils(iq) 1009 endif !if (nqfils(iq).gt.0) then 1010 ! end CRisi 1011 860 1012 DO l=1,llm 861 1013 DO ij=1,ip1jmp1 862 newmasse=masse(ij,l )+w(ij,l+1)-w(ij,l)863 q(ij,l )=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))1014 newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l) 1015 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l)) 864 1016 & /newmasse 865 masse(ij,l)=newmasse 866 ENDDO 867 ENDDO 868 1017 masse(ij,l,iq)=newmasse 1018 ENDDO 1019 ENDDO 1020 1021 ! retablir les fils en rapport de melange par rapport a l'air: 1022 if (nqfils(iq).gt.0) then 1023 do ifils=1,nqdesc(iq) 1024 iq2=iqfils(ifils,iq) 1025 DO l=1,llm 1026 DO ij=1,ip1jmp1 1027 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 1028 enddo 1029 enddo 1030 enddo !do ifils=1,nqdesc(iq) 1031 endif !if (nqfils(iq).gt.0) then 1032 !write(*,*) 'vlsplt 1032' 869 1033 870 1034 RETURN -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d/vlspltqs.F ¶
r1422 r1508 1 ! 2 ! $Header$3 ! 1 c 2 c $Id: vlspltqs.F 2286 2015-05-20 13:27:07Z emillour $ 3 c 4 4 SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt, 5 , p,pk,teta ) 5 , p,pk,teta,iq ) 6 USE infotrac, ONLY: nqtot,nqdesc,iqfils 6 7 c 7 8 c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron … … 33 34 REAL masse(ip1jmp1,llm),pente_max 34 35 REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm) 35 REAL q(ip1jmp1,llm )36 REAL q(ip1jmp1,llm,nqtot) 36 37 REAL w(ip1jmp1,llm),pdt 37 38 REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm) 39 INTEGER iq ! CRisi 38 40 c 39 41 c Local … … 41 43 c 42 44 INTEGER i,ij,l,j,ii 45 INTEGER ifils,iq2 ! CRisi 43 46 c 44 47 REAL qsat(ip1jmp1,llm) 45 REAL zm(ip1jmp1,llm )48 REAL zm(ip1jmp1,llm,nqtot) 46 49 REAL mu(ip1jmp1,llm) 47 50 REAL mv(ip1jm,llm) 48 51 REAL mw(ip1jmp1,llm+1) 49 REAL zq(ip1jmp1,llm )52 REAL zq(ip1jmp1,llm,nqtot) 50 53 REAL temps1,temps2,temps3 51 54 REAL zzpbar, zzw … … 63 66 REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play 64 67 REAL ptarg,pdelarg,foeew,zdelta 65 ! ADAPTATION GCM POUR CP(T)66 68 REAL tempe(ip1jmp1,llm) 67 69 … … 115 117 ENDDO 116 118 117 CALL SCOPY(ijp1llm,q,1,zq,1) 118 CALL SCOPY(ijp1llm,masse,1,zm,1) 119 CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1) 120 CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1) 121 if (nqdesc(iq).gt.0) then 122 do ifils=1,nqdesc(iq) 123 iq2=iqfils(ifils,iq) 124 CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1) 125 enddo 126 endif !if (nqfils(iq).gt.0) then 119 127 120 128 c call minmaxq(zq,qmin,qmax,'avant vlxqs ') 121 call vlxqs(zq,pente_max,zm,mu,qsat) 122 129 call vlxqs(zq,pente_max,zm,mu,qsat,iq) 123 130 124 131 c call minmaxq(zq,qmin,qmax,'avant vlyqs ') 125 132 126 call vlyqs(zq,pente_max,zm,mv,qsat) 127 133 call vlyqs(zq,pente_max,zm,mv,qsat,iq) 128 134 129 135 c call minmaxq(zq,qmin,qmax,'avant vlz ') 130 136 131 call vlz(zq,pente_max,zm,mw) 132 137 call vlz(zq,pente_max,zm,mw,iq) 133 138 134 139 c call minmaxq(zq,qmin,qmax,'avant vlyqs ') 135 140 c call minmaxq(zm,qmin,qmax,'M avant vlyqs ') 136 141 137 call vlyqs(zq,pente_max,zm,mv,qsat) 138 142 call vlyqs(zq,pente_max,zm,mv,qsat,iq) 139 143 140 144 c call minmaxq(zq,qmin,qmax,'avant vlxqs ') 141 145 c call minmaxq(zm,qmin,qmax,'M avant vlxqs ') 142 146 143 call vlxqs(zq,pente_max,zm,mu,qsat )147 call vlxqs(zq,pente_max,zm,mu,qsat,iq) 144 148 145 149 c call minmaxq(zq,qmin,qmax,'apres vlxqs ') … … 149 153 DO l=1,llm 150 154 DO ij=1,ip1jmp1 151 q(ij,l )=zq(ij,l)155 q(ij,l,iq)=zq(ij,l,iq) 152 156 ENDDO 153 157 DO ij=1,ip1jm+1,iip1 154 q(ij+iim,l)=q(ij,l) 155 ENDDO 156 ENDDO 158 q(ij+iim,l,iq)=q(ij,l,iq) 159 ENDDO 160 ENDDO 161 ! CRisi: aussi pour les fils 162 if (nqdesc(iq).gt.0) then 163 do ifils=1,nqdesc(iq) 164 iq2=iqfils(ifils,iq) 165 DO l=1,llm 166 DO ij=1,ip1jmp1 167 q(ij,l,iq2)=zq(ij,l,iq2) 168 ENDDO 169 DO ij=1,ip1jm+1,iip1 170 q(ij+iim,l,iq2)=q(ij,l,iq2) 171 ENDDO 172 ENDDO 173 enddo !do ifils=1,nqdesc(iq) 174 endif ! if (nqfils(iq).gt.0) then 175 !write(*,*) 'vlspltqs 183: fin de la routine' 157 176 158 177 RETURN 159 178 END 160 SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat) 179 SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq) 180 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 181 161 182 c 162 183 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 175 196 c Arguments: 176 197 c ---------- 177 REAL masse(ip1jmp1,llm ),pente_max198 REAL masse(ip1jmp1,llm,nqtot),pente_max 178 199 REAL u_m( ip1jmp1,llm ) 179 REAL q(ip1jmp1,llm )200 REAL q(ip1jmp1,llm,nqtot) 180 201 REAL qsat(ip1jmp1,llm) 202 INTEGER iq ! CRisi 181 203 c 182 204 c Local … … 191 213 REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm) 192 214 REAL u_mq(ip1jmp1,llm) 215 216 ! CRisi 217 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) 218 INTEGER ifils,iq2 ! CRisi 193 219 194 220 Logical first,testcpu … … 223 249 DO l = 1, llm 224 250 DO ij=iip2,ip1jm-1 225 dxqu(ij)=q(ij+1,l )-q(ij,l)251 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 226 252 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0' 227 c sigu(ij)=u_m(ij,l)/masse(ij,l )253 c sigu(ij)=u_m(ij,l)/masse(ij,l,iq) 228 254 ENDDO 229 255 DO ij=iip1+iip1,ip1jm,iip1 … … 277 303 DO l = 1, llm 278 304 DO ij=iip2,ip1jm-1 279 dxqu(ij)=q(ij+1,l )-q(ij,l)305 dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq) 280 306 ENDDO 281 307 DO ij=iip1+iip1,ip1jm,iip1 … … 319 345 DO l=1,llm 320 346 DO ij=iip2,ip1jm-1 321 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l ),322 , 1.+u_m(ij,l)/masse(ij+1,l ),347 zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq), 348 , 1.+u_m(ij,l)/masse(ij+1,l,iq), 323 349 , u_m(ij,l)) 324 350 zdum(ij,l)=0.5*zdum(ij,l) 325 351 u_mq(ij,l)=cvmgp( 326 , q(ij,l )+zdum(ij,l)*dxq(ij,l),327 , q(ij+1,l )-zdum(ij,l)*dxq(ij+1,l),352 , q(ij,l,iq)+zdum(ij,l)*dxq(ij,l), 353 , q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l), 328 354 , u_m(ij,l)) 329 355 u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l) … … 337 363 DO ij=iip2,ip1jm-1 338 364 IF (u_m(ij,l).gt.0.) THEN 339 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l )365 zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq) 340 366 u_mq(ij,l)=u_m(ij,l)* 341 $ min(q(ij,l )+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))367 $ min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l)) 342 368 ELSE 343 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l )369 zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq) 344 370 u_mq(ij,l)=u_m(ij,l)* 345 $ min(q(ij+1,l )-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))371 $ min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l)) 346 372 ENDIF 347 373 ENDDO … … 412 438 i=ijq-(j-1)*iip1 413 439 c accumulation pour les mailles completements advectees 414 do while(zu_m.gt.masse(ijq,l)) 415 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l) 416 zu_m=zu_m-masse(ijq,l) 440 do while(zu_m.gt.masse(ijq,l,iq)) 441 u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq) 442 & *masse(ijq,l,iq) 443 zu_m=zu_m-masse(ijq,l,iq) 417 444 i=mod(i-2+iim,iim)+1 418 445 ijq=(j-1)*iip1+i … … 420 447 c ajout de la maille non completement advectee 421 448 u_mq(ij,l)=u_mq(ij,l)+zu_m* 422 & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l)) 449 & (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq)) 450 & *dxq(ijq,l)) 423 451 ELSE 424 452 ijq=ij+1 425 453 i=ijq-(j-1)*iip1 426 454 c accumulation pour les mailles completements advectees 427 do while(-zu_m.gt.masse(ijq,l)) 428 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l) 429 zu_m=zu_m+masse(ijq,l) 455 do while(-zu_m.gt.masse(ijq,l,iq)) 456 u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) 457 & *masse(ijq,l,iq) 458 zu_m=zu_m+masse(ijq,l,iq) 430 459 i=mod(i,iim)+1 431 460 ijq=(j-1)*iip1+i 432 461 ENDDO 433 462 c ajout de la maille non completement advectee 434 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l )-435 & 0.5*(1.+zu_m/masse(ijq,l ))*dxq(ijq,l))463 u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- 464 & 0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l)) 436 465 ENDIF 437 466 ENDDO … … 450 479 ENDDO 451 480 481 ! CRisi: appel récursif de l'advection sur les fils. 482 ! Il faut faire ça avant d'avoir mis à jour q et masse 483 !write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq) 484 485 if (nqfils(iq).gt.0) then 486 do ifils=1,nqdesc(iq) 487 iq2=iqfils(ifils,iq) 488 DO l=1,llm 489 DO ij=iip2,ip1jm 490 ! On a besoin de q et masse seulement entre iip2 et ip1jm 491 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 492 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 493 enddo 494 enddo 495 enddo !do ifils=1,nqdesc(iq) 496 do ifils=1,nqfils(iq) 497 iq2=iqfils(ifils,iq) 498 call vlx(Ratio,pente_max,masseq,u_mq,iq2) 499 enddo !do ifils=1,nqfils(iq) 500 endif !if (nqfils(iq).gt.0) then 501 ! end CRisi 452 502 453 503 c calcul des tendances … … 455 505 DO l=1,llm 456 506 DO ij=iip2+1,ip1jm 457 new_m=masse(ij,l )+u_m(ij-1,l)-u_m(ij,l)458 q(ij,l )=(q(ij,l)*masse(ij,l)+507 new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l) 508 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ 459 509 & u_mq(ij-1,l)-u_mq(ij,l)) 460 510 & /new_m 461 masse(ij,l )=new_m511 masse(ij,l,iq)=new_m 462 512 ENDDO 463 513 c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous) 464 514 DO ij=iip1+iip1,ip1jm,iip1 465 q(ij-iim,l)=q(ij,l) 466 masse(ij-iim,l)=masse(ij,l) 467 ENDDO 468 ENDDO 515 q(ij-iim,l,iq)=q(ij,l,iq) 516 masse(ij-iim,l,iq)=masse(ij,l,iq) 517 ENDDO 518 ENDDO 519 520 ! retablir les fils en rapport de melange par rapport a l'air: 521 ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio 522 ! puis on boucle en longitude 523 if (nqdesc(iq).gt.0) then 524 do ifils=1,nqdesc(iq) 525 iq2=iqfils(ifils,iq) 526 DO l=1,llm 527 DO ij=iip2+1,ip1jm 528 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 529 enddo 530 DO ij=iip1+iip1,ip1jm,iip1 531 q(ij-iim,l,iq2)=q(ij,l,iq2) 532 enddo ! DO ij=ijb+iip1-1,ije,iip1 533 enddo !DO l=1,llm 534 enddo !do ifils=1,nqdesc(iq) 535 endif !if (nqfils(iq).gt.0) then 469 536 470 537 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1) … … 474 541 RETURN 475 542 END 476 SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat) 543 SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq) 544 USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 545 USE comconst_mod, ONLY: pi 477 546 c 478 547 c Auteurs: P.Le Van, F.Hourdin, F.Forget … … 486 555 c 487 556 c -------------------------------------------------------------------- 488 USE comconst_mod, ONLY: pi489 490 557 IMPLICIT NONE 491 558 c … … 497 564 c Arguments: 498 565 c ---------- 499 REAL masse(ip1jmp1,llm ),pente_max566 REAL masse(ip1jmp1,llm,nqtot),pente_max 500 567 REAL masse_adv_v( ip1jm,llm) 501 REAL q(ip1jmp1,llm )568 REAL q(ip1jmp1,llm,nqtot) 502 569 REAL qsat(ip1jmp1,llm) 570 INTEGER iq ! CRisi 503 571 c 504 572 c Local … … 524 592 SAVE sinlon,coslon,sinlondlon,coslondlon 525 593 SAVE airej2,airejjm 594 595 REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi 596 INTEGER ifils,iq2 ! CRisi 526 597 c 527 598 c … … 562 633 563 634 DO i = 1, iim 564 airescb(i) = aire(i+ iip1) * q(i+ iip1,l )565 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l )635 airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq) 636 airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq) 566 637 ENDDO 567 638 qpns = SSUM( iim, airescb ,1 ) / airej2 … … 571 642 572 643 DO ij=1,ip1jm 573 dyqv(ij)=q(ij,l )-q(ij+iip1,l)644 dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq) 574 645 adyqv(ij)=abs(dyqv(ij)) 575 646 ENDDO … … 586 657 587 658 DO ij=1,iip1 588 dyq(ij,l)=qpns-q(ij+iip1,l )589 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l )-qpsn659 dyq(ij,l)=qpns-q(ij+iip1,l,iq) 660 dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn 590 661 ENDDO 591 662 … … 705 776 DO ij=1,ip1jm 706 777 IF( masse_adv_v(ij,l).GT.0. ) THEN 707 qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) + 708 , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))) 778 qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq ) + 779 , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l) 780 , /masse(ij+iip1,l,iq))) 709 781 ELSE 710 qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l ) - dyq(ij,l) *711 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l )) )782 qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) * 783 , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) ) 712 784 ENDIF 713 785 qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l) … … 716 788 717 789 790 ! CRisi: appel récursif de l'advection sur les fils. 791 ! Il faut faire ça avant d'avoir mis à jour q et masse 792 !write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq) 793 794 if (nqfils(iq).gt.0) then 795 do ifils=1,nqdesc(iq) 796 iq2=iqfils(ifils,iq) 797 DO l=1,llm 798 DO ij=1,ip1jmp1 799 masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq) 800 Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq) 801 enddo 802 enddo 803 enddo !do ifils=1,nqdesc(iq) 804 805 do ifils=1,nqfils(iq) 806 iq2=iqfils(ifils,iq) 807 !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2 808 call vly(Ratio,pente_max,masseq,qbyv,iq2) 809 enddo !do ifils=1,nqfils(iq) 810 endif !if (nqfils(iq).gt.0) then 811 718 812 DO l=1,llm 719 813 DO ij=iip2,ip1jm 720 newmasse=masse(ij,l )814 newmasse=masse(ij,l,iq) 721 815 & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l) 722 q(ij,l )=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))723 & /newmasse724 masse(ij,l )=newmasse816 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l) 817 & -qbyv(ij-iip1,l))/newmasse 818 masse(ij,l,iq)=newmasse 725 819 ENDDO 726 820 c.-. ancienne version … … 728 822 convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln 729 823 DO ij = 1,iip1 730 newmasse=masse(ij,l )+convmpn*aire(ij)731 q(ij,l )=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/824 newmasse=masse(ij,l,iq)+convmpn*aire(ij) 825 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/ 732 826 & newmasse 733 masse(ij,l )=newmasse827 masse(ij,l,iq)=newmasse 734 828 ENDDO 735 829 convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols 736 830 convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols 737 831 DO ij = ip1jm+1,ip1jmp1 738 newmasse=masse(ij,l )+convmps*aire(ij)739 q(ij,l )=(q(ij,l)*masse(ij,l)+convps*aire(ij))/832 newmasse=masse(ij,l,iq)+convmps*aire(ij) 833 q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/ 740 834 & newmasse 741 masse(ij,l )=newmasse835 masse(ij,l,iq)=newmasse 742 836 ENDDO 743 837 c.-. fin ancienne version … … 752 846 c DO ij = 1,iip1 753 847 c q(ij,l)=newq 754 c masse(ij,l )=newmasse*aire(ij)848 c masse(ij,l,iq)=newmasse*aire(ij) 755 849 c ENDDO 756 850 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1) … … 762 856 c DO ij = ip1jm+1,ip1jmp1 763 857 c q(ij,l)=newq 764 c masse(ij,l )=newmasse*aire(ij)858 c masse(ij,l,iq)=newmasse*aire(ij) 765 859 c ENDDO 766 860 c._. fin nouvelle version 767 861 ENDDO 768 862 863 !write(*,*) 'vly 866' 864 865 ! retablir les fils en rapport de melange par rapport a l'air: 866 if (nqdesc(iq).gt.0) then 867 do ifils=1,nqdesc(iq) 868 iq2=iqfils(ifils,iq) 869 DO l=1,llm 870 DO ij=1,ip1jmp1 871 q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2) 872 enddo 873 enddo 874 enddo !do ifils=1,nqdesc(iq) 875 endif !if (nqfils(iq).gt.0) then 876 !write(*,*) 'vly 879' 877 769 878 RETURN 770 879 END -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d_common/covcont.F90 ¶
r1506 r1508 1 SUBROUTINE covcont (klevel,ucov, vcov, ucont, vcont ) 1 2 ! 2 ! $Header$ 3 ! 4 SUBROUTINE covcont (klevel,ucov, vcov, ucont, vcont ) 5 IMPLICIT NONE 3 !------------------------------------------------------------------------------- 4 ! Author: P. Le Van 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute contravariant components from covariant components. 7 !------------------------------------------------------------------------------- 8 IMPLICIT NONE 9 include "dimensions.h" 10 include "paramet.h" 11 include "comgeom.h" 12 !=============================================================================== 13 ! Arguments: 14 INTEGER, INTENT(IN) :: klevel !--- VERTICAL LEVELS NUMBER 15 REAL, INTENT(IN) :: ucov ( ip1jmp1,klevel ) !--- U COVARIANT WIND 16 REAL, INTENT(IN) :: vcov ( ip1jm ,klevel ) !--- V COVARIANT WIND 17 REAL, INTENT(OUT) :: ucont( ip1jmp1,klevel ) !--- U CONTRAVAR WIND 18 REAL, INTENT(OUT) :: vcont( ip1jm ,klevel ) !--- V CONTRAVAR WIND 19 !=============================================================================== 20 ! Local variables: 21 INTEGER :: l 22 !=============================================================================== 23 DO l=1,klevel 24 ucont(iip2:ip1jm,l)=ucov(iip2:ip1jm,l) * unscu2(iip2:ip1jm) 25 vcont( 1:ip1jm,l)=vcov( 1:ip1jm,l) * unscv2( 1:ip1jm) 26 END DO 6 27 7 c======================================================================= 8 c 9 c Auteur: P. Le Van 10 c ------- 11 c 12 c Objet: 13 c ------ 14 c 15 c ********************************************************************* 16 c calcul des compos. contravariantes a partir des comp.covariantes 17 c ******************************************************************** 18 c 19 c======================================================================= 28 END SUBROUTINE covcont 20 29 21 #include "dimensions.h"22 #include "paramet.h"23 #include "comgeom.h"24 25 INTEGER klevel26 REAL ucov( ip1jmp1,klevel ), vcov( ip1jm,klevel )27 REAL ucont( ip1jmp1,klevel ), vcont( ip1jm,klevel )28 INTEGER l,ij29 30 31 DO 10 l = 1,klevel32 33 DO 2 ij = iip2, ip1jm34 ucont( ij,l ) = ucov( ij,l ) * unscu2( ij )35 2 CONTINUE36 37 DO 4 ij = 1,ip1jm38 vcont( ij,l ) = vcov( ij,l ) * unscv2( ij )39 4 CONTINUE40 41 10 CONTINUE42 RETURN43 END -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d_common/dynetat0.F90 ¶
r1506 r1508 2 2 ! $Id $ 3 3 ! 4 SUBROUTINE dynetat0(fichnom,vcov,ucov, 5 . teta,q,masse,ps,phis,time0) 6 7 USE infotrac, only: tname, nqtot 8 use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, 9 & nf90_get_var, nf90_inq_varid, nf90_inq_dimid, 10 & nf90_inquire_dimension,nf90_close 4 SUBROUTINE dynetat0(fichnom,vcov,ucov,teta,q,masse,ps,phis,time0) 5 6 USE infotrac, only: tname, nqtot, zone_num, iso_indnum,& 7 iso_num, phase_num, alpha_ideal, iqiso, & 8 ok_isotopes, iqpere, tnat 9 use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, & 10 nf90_get_var, nf90_inq_varid, nf90_inq_dimid, & 11 nf90_inquire_dimension,nf90_close 11 12 12 13 use control_mod, only : planet_type, timestart 13 14 USE comvert_mod, ONLY: pa,preff 14 USE comconst_mod, ONLY: im,jm,lllm,daysec,dtvr, 15 .rad,omeg,g,cpp,kappa,pi15 USE comconst_mod, ONLY: im,jm,lllm,daysec,dtvr, & 16 rad,omeg,g,cpp,kappa,pi 16 17 USE logic_mod, ONLY: fxyhypb,ysinus 17 18 USE serre_mod, ONLY: clon,clat,grossismx,grossismy 18 USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn, 19 .start_time,day_ini,hour_ini19 USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn, & 20 start_time,day_ini,hour_ini 20 21 USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 21 22 22 23 IMPLICIT NONE 23 24 24 c======================================================================= 25 c 26 c Auteur: P. Le Van / L.Fairhead 27 c ------- 28 c 29 c objet: 30 c ------ 31 c 32 c Lecture de l'etat initial 33 c 34 c======================================================================= 35 c----------------------------------------------------------------------- 36 c Declarations: 37 c ------------- 38 39 #include "dimensions.h" 40 #include "paramet.h" 41 #include "comgeom2.h" 42 #include "netcdf.inc" 43 #include "iniprint.h" 44 45 c Arguments: 46 c ---------- 47 48 CHARACTER(len=*),INTENT(IN) :: fichnom 49 REAL,INTENT(OUT) :: vcov(iip1,jjm,llm) 50 REAL,INTENT(OUT) :: ucov(iip1,jjp1,llm) 51 REAL,INTENT(OUT) :: teta(iip1,jjp1,llm) 52 REAL,INTENT(OUT) :: q(iip1,jjp1,llm,nqtot) 53 REAL,INTENT(OUT) :: masse(iip1,jjp1,llm) 54 REAL,INTENT(OUT) :: ps(iip1,jjp1) 55 REAL,INTENT(OUT) :: phis(iip1,jjp1) 56 REAL,INTENT(OUT) :: time0 57 58 c Variables 59 c 60 INTEGER length,iq 61 PARAMETER (length = 100) 62 REAL tab_cntrl(length) ! tableau des parametres du run 63 INTEGER ierr, nid, nvarid 64 65 character(len=12) :: start_file_type="earth" ! default start file type 66 INTEGER idecal 67 68 69 REAL,ALLOCATABLE :: time(:) ! times stored in start 70 INTEGER timelen ! number of times stored in the file 71 INTEGER indextime ! index of selected time 72 !REAL hour_ini ! fraction of day of stored date. Equivalent of day_ini, but 0=<hour_ini<1 73 74 INTEGER edges(4),corner(4) 75 integer :: i 76 77 c----------------------------------------------------------------------- 78 79 c Ouverture NetCDF du fichier etat initial 80 81 ierr=nf90_open(fichnom,NF90_NOWRITE,nid) 82 IF (ierr.NE.nf90_noerr) THEN 83 write(lunout,*)'dynetat0: Pb d''ouverture du fichier start.nc' 84 write(lunout,*)trim(nf90_strerror(ierr)) 85 CALL ABORT_gcm("dynetat0", "", 1) 86 ENDIF 87 88 c 89 ierr = nf90_inq_varid (nid, "controle", nvarid) 90 IF (ierr .NE. nf90_noerr) THEN 91 write(lunout,*)"dynetat0: Le champ <controle> est absent" 92 write(lunout,*)trim(nf90_strerror(ierr)) 93 CALL ABORT_gcm("dynetat0", "", 1) 94 ENDIF 95 ierr = nf90_get_var(nid, nvarid, tab_cntrl) 96 IF (ierr .NE. nf90_noerr) THEN 97 write(lunout,*)"dynetat0: Lecture echoue pour <controle>" 98 write(lunout,*)trim(nf90_strerror(ierr)) 99 CALL ABORT_gcm("dynetat0", "", 1) 100 ENDIF 25 !======================================================================= 26 ! 27 ! Read initial confitions file 28 ! 29 !======================================================================= 30 31 include "dimensions.h" 32 include "paramet.h" 33 include "comgeom2.h" 34 include "iniprint.h" 35 36 !=============================================================================== 37 ! Arguments: 38 CHARACTER(LEN=*), INTENT(IN) :: fichnom !--- FILE NAME 39 REAL, INTENT(OUT) :: vcov(iip1,jjm, llm) !--- V COVARIANT WIND 40 REAL, INTENT(OUT) :: ucov(iip1,jjp1,llm) !--- U COVARIANT WIND 41 REAL, INTENT(OUT) :: teta(iip1,jjp1,llm) !--- POTENTIAL TEMP. 42 REAL, INTENT(OUT) :: q(iip1,jjp1,llm,nqtot) !--- TRACERS 43 REAL, INTENT(OUT) :: masse(iip1,jjp1,llm) !--- MASS PER CELL 44 REAL, INTENT(OUT) :: ps(iip1,jjp1) !--- GROUND PRESSURE 45 REAL, INTENT(OUT) :: phis(iip1,jjp1) !--- GEOPOTENTIAL 46 REAL,INTENT(OUT) :: time0 47 !=============================================================================== 48 ! Local Variables 49 CHARACTER(LEN=256) :: msg, var, modname 50 INTEGER,PARAMETER :: length=100 51 INTEGER :: iq, fID, vID, idecal 52 REAL :: tab_cntrl(length) ! array containing run parameters 53 INTEGER :: ierr 54 CHARACTER(len=12) :: start_file_type="earth" ! default start file type 55 56 REAL,ALLOCATABLE :: time(:) ! times stored in start 57 INTEGER :: timelen ! number of times stored in the file 58 INTEGER :: indextime ! index of selected time 59 !REAL hour_ini ! fraction of day of stored date. Equivalent of day_ini, but 0=<hour_ini<1 60 61 INTEGER :: edges(4),corner(4) 62 INTEGER :: i 63 64 !----------------------------------------------------------------------- 65 modname="dynetat0" 66 67 ! Open initial state NetCDF file 68 var=fichnom 69 CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var) 70 ! 71 CALL get_var1("controle",tab_cntrl) 101 72 102 73 !!! AS: idecal is a hack to be able to read planeto starts... … … 137 108 pa = tab_cntrl(idecal+13) 138 109 preff = tab_cntrl(idecal+14) 139 c 110 ! 140 111 clon = tab_cntrl(idecal+15) 141 112 clat = tab_cntrl(idecal+16) 142 113 grossismx = tab_cntrl(idecal+17) 143 114 grossismy = tab_cntrl(idecal+18) 144 c 115 ! 145 116 IF ( tab_cntrl(idecal+19).EQ.1. ) THEN 146 fxyhypb = . TRUE.147 cdzoomx = tab_cntrl(25)148 cdzoomy = tab_cntrl(26)149 ctaux = tab_cntrl(28)150 ctauy = tab_cntrl(29)117 fxyhypb = .TRUE. 118 ! dzoomx = tab_cntrl(25) 119 ! dzoomy = tab_cntrl(26) 120 ! taux = tab_cntrl(28) 121 ! tauy = tab_cntrl(29) 151 122 ELSE 152 fxyhypb = . FALSE.153 ysinus = . FALSE.154 IF( tab_cntrl(idecal+22).EQ.1. ) ysinus = . 123 fxyhypb = .FALSE. 124 ysinus = .FALSE. 125 IF( tab_cntrl(idecal+22).EQ.1. ) ysinus = .TRUE. 155 126 ENDIF 156 127 … … 170 141 start_time=0 171 142 endif 172 c ................................................................. 173 c 174 c 175 write(lunout,*)'dynetat0: rad,omeg,g,cpp,kappa ', 176 & rad,omeg,g,cpp,kappa 177 178 IF( im.ne.iim ) THEN 179 PRINT 1,im,iim 180 STOP 181 ELSE IF( jm.ne.jjm ) THEN 182 PRINT 2,jm,jjm 183 STOP 184 ELSE IF( lllm.ne.llm ) THEN 185 PRINT 3,lllm,llm 186 STOP 187 ENDIF 188 189 ierr=nf90_inq_varid(nid, "rlonu", nvarid) 190 IF (ierr .NE. nf90_noerr) THEN 191 write(lunout,*)"dynetat0: Le champ <rlonu> est absent" 192 write(lunout,*)trim(nf90_strerror(ierr)) 193 CALL ABORT_gcm("dynetat0", "", 1) 194 ENDIF 195 ierr = nf90_get_var(nid, nvarid, rlonu) 196 IF (ierr .NE. nf90_noerr) THEN 197 write(lunout,*)"dynetat0: Lecture echouee pour <rlonu>" 198 write(lunout,*)trim(nf90_strerror(ierr)) 199 CALL ABORT_gcm("dynetat0", "", 1) 200 ENDIF 201 202 ierr = nf90_inq_varid (nid, "rlatu", nvarid) 203 IF (ierr .NE. nf90_noerr) THEN 204 write(lunout,*)"dynetat0: Le champ <rlatu> est absent" 205 write(lunout,*)trim(nf90_strerror(ierr)) 206 CALL ABORT_gcm("dynetat0", "", 1) 207 ENDIF 208 ierr = nf90_get_var(nid, nvarid, rlatu) 209 IF (ierr .NE. nf90_noerr) THEN 210 write(lunout,*)"dynetat0: Lecture echouee pour <rlatu>" 211 write(lunout,*)trim(nf90_strerror(ierr)) 212 CALL ABORT_gcm("dynetat0", "", 1) 213 ENDIF 214 215 ierr = nf90_inq_varid (nid, "rlonv", nvarid) 216 IF (ierr .NE. nf90_noerr) THEN 217 write(lunout,*)"dynetat0: Le champ <rlonv> est absent" 218 write(lunout,*)trim(nf90_strerror(ierr)) 219 CALL ABORT_gcm("dynetat0", "", 1) 220 ENDIF 221 ierr = nf90_get_var(nid, nvarid, rlonv) 222 IF (ierr .NE. nf90_noerr) THEN 223 write(lunout,*)"dynetat0: Lecture echouee pour <rlonv>" 224 write(lunout,*)trim(nf90_strerror(ierr)) 225 CALL ABORT_gcm("dynetat0", "", 1) 226 ENDIF 227 228 ierr = nf90_inq_varid (nid, "rlatv", nvarid) 229 IF (ierr .NE. nf90_noerr) THEN 230 write(lunout,*)"dynetat0: Le champ <rlatv> est absent" 231 write(lunout,*)trim(nf90_strerror(ierr)) 232 CALL ABORT_gcm("dynetat0", "", 1) 233 ENDIF 234 ierr = nf90_get_var(nid, nvarid, rlatv) 235 IF (ierr .NE. nf90_noerr) THEN 236 write(lunout,*)"dynetat0: Lecture echouee pour rlatv" 237 write(lunout,*)trim(nf90_strerror(ierr)) 238 CALL ABORT_gcm("dynetat0", "", 1) 239 ENDIF 240 241 ierr = nf90_inq_varid (nid, "cu", nvarid) 242 IF (ierr .NE. nf90_noerr) THEN 243 write(lunout,*)"dynetat0: Le champ <cu> est absent" 244 write(lunout,*)trim(nf90_strerror(ierr)) 245 CALL ABORT_gcm("dynetat0", "", 1) 246 ENDIF 247 ierr = nf90_get_var(nid, nvarid, cu) 248 IF (ierr .NE. nf90_noerr) THEN 249 write(lunout,*)"dynetat0: Lecture echouee pour <cu>" 250 write(lunout,*)trim(nf90_strerror(ierr)) 251 CALL ABORT_gcm("dynetat0", "", 1) 252 ENDIF 253 254 ierr = nf90_inq_varid (nid, "cv", nvarid) 255 IF (ierr .NE. nf90_noerr) THEN 256 write(lunout,*)"dynetat0: Le champ <cv> est absent" 257 write(lunout,*)trim(nf90_strerror(ierr)) 258 CALL ABORT_gcm("dynetat0", "", 1) 259 ENDIF 260 ierr = nf90_get_var(nid, nvarid, cv) 261 IF (ierr .NE. nf90_noerr) THEN 262 write(lunout,*)"dynetat0: Lecture echouee pour <cv>" 263 write(lunout,*)trim(nf90_strerror(ierr)) 264 CALL ABORT_gcm("dynetat0", "", 1) 265 ENDIF 266 267 ierr = nf90_inq_varid (nid, "aire", nvarid) 268 IF (ierr .NE. nf90_noerr) THEN 269 write(lunout,*)"dynetat0: Le champ <aire> est absent" 270 write(lunout,*)trim(nf90_strerror(ierr)) 271 CALL ABORT_gcm("dynetat0", "", 1) 272 ENDIF 273 ierr = nf90_get_var(nid, nvarid, aire) 274 IF (ierr .NE. nf90_noerr) THEN 275 write(lunout,*)"dynetat0: Lecture echouee pour <aire>" 276 write(lunout,*)trim(nf90_strerror(ierr)) 277 CALL ABORT_gcm("dynetat0", "", 1) 278 ENDIF 279 280 ierr = nf90_inq_varid (nid, "phisinit", nvarid) 281 IF (ierr .NE. nf90_noerr) THEN 282 write(lunout,*)"dynetat0: Le champ <phisinit> est absent" 283 write(lunout,*)trim(nf90_strerror(ierr)) 284 CALL ABORT_gcm("dynetat0", "", 1) 285 ENDIF 286 ierr = nf90_get_var(nid, nvarid, phis) 287 IF (ierr .NE. nf90_noerr) THEN 288 write(lunout,*)"dynetat0: Lecture echouee pour <phisinit>" 289 write(lunout,*)trim(nf90_strerror(ierr)) 290 CALL ABORT_gcm("dynetat0", "", 1) 291 ENDIF 143 ! ................................................................. 144 ! 145 ! 146 WRITE(lunout,*)trim(modname)//': rad,omeg,g,cpp,kappa ', & 147 rad,omeg,g,cpp,kappa 148 149 CALL check_dim(im,iim,'im','iim') 150 CALL check_dim(jm,jjm,'jm','jjm') 151 CALL check_dim(llm,lllm,'llm','lllm') 152 153 CALL get_var1("rlonu",rlonu) 154 CALL get_var1("rlatu",rlatu) 155 CALL get_var1("rlonv",rlonv) 156 CALL get_var1("rlatv",rlatv) 157 158 CALL get_var2("cu" ,cu) 159 CALL get_var2("cv" ,cv) 160 161 CALL get_var2("aire" ,aire) 162 CALL get_var2("phisinit",phis) 292 163 293 164 ! read time axis 294 ierr = nf90_inq_varid ( nid, "temps", nvarid)165 ierr = nf90_inq_varid (fID, "temps", vID) 295 166 IF (ierr .NE. nf90_noerr) THEN 296 167 write(lunout,*)"dynetat0: Le champ <temps> est absent" 297 168 write(lunout,*)"dynetat0: J essaie <Time>" 298 ierr = nf90_inq_varid ( nid, "Time", nvarid)169 ierr = nf90_inq_varid (fID, "Time", vID) 299 170 IF (ierr .NE. nf90_noerr) THEN 300 171 write(lunout,*)"dynetat0: Le champ <Time> est absent" … … 303 174 ENDIF 304 175 ! Get the length of the "Time" dimension 305 ierr = nf90_inq_dimid( nid,"Time",nvarid)306 ierr = nf90_inquire_dimension( nid,nvarid,len=timelen)176 ierr = nf90_inq_dimid(fID,"Time",vID) 177 ierr = nf90_inquire_dimension(fID,vID,len=timelen) 307 178 allocate(time(timelen)) 308 179 ! Then look for the "Time" variable 309 ierr =nf90_inq_varid( nid,"Time",nvarid)310 ierr = nf90_get_var( nid, nvarid, time)180 ierr =nf90_inq_varid(fID,"Time",vID) 181 ierr = nf90_get_var(fID, vID, time) 311 182 IF (ierr .NE. nf90_noerr) THEN 312 183 write(lunout,*)"dynetat0: Lecture echouee <Time>" … … 316 187 ELSE 317 188 ! Get the length of the "temps" dimension 318 ierr = nf90_inq_dimid( nid,"temps",nvarid)319 ierr = nf90_inquire_dimension( nid,nvarid,len=timelen)189 ierr = nf90_inq_dimid(fID,"temps",vID) 190 ierr = nf90_inquire_dimension(fID,vID,len=timelen) 320 191 allocate(time(timelen)) 321 192 ! Then look for the "temps" variable 322 ierr = nf90_inq_varid ( nid, "temps", nvarid)323 ierr = nf90_get_var( nid, nvarid, time)193 ierr = nf90_inq_varid (fID, "temps", vID) 194 ierr = nf90_get_var(fID, vID, time) 324 195 IF (ierr .NE. nf90_noerr) THEN 325 196 write(lunout,*)"dynetat0: Lecture echouee <temps>" … … 341 212 ENDDO 342 213 IF (indextime .eq. 0) THEN 343 write(lunout,*)"Time", timestart," is not in " 344 &//trim(fichnom)//"!!"214 write(lunout,*)"Time", timestart," is not in " & 215 //trim(fichnom)//"!!" 345 216 write(lunout,*)"Stored times are:" 346 217 DO i=1,timelen … … 362 233 endif 363 234 364 PRINT*, "dynetat0: Selected time ",time(indextime), 365 ." at index ",indextime235 PRINT*, "dynetat0: Selected time ",time(indextime), & 236 " at index ",indextime 366 237 367 238 DEALLOCATE(time) 368 239 369 240 ! read vcov 370 corner(1)=1 371 corner(2)=1 372 corner(3)=1 373 corner(4)=indextime 374 edges(1)=iip1 375 edges(2)=jjm 376 edges(3)=llm 377 edges(4)=1 378 ierr=nf90_inq_varid(nid,"vcov",nvarid) 379 IF (ierr .NE. nf90_noerr) THEN 380 write(lunout,*)"dynetat0: Le champ <vcov> est absent" 381 write(lunout,*)trim(nf90_strerror(ierr)) 382 CALL ABORT_gcm("dynetat0", "", 1) 383 ENDIF 384 ierr=nf90_get_var(nid,nvarid,vcov,corner,edges) 385 IF (ierr .NE. nf90_noerr) THEN 386 write(lunout,*)"dynetat0: Lecture echouee pour <vcov>" 387 write(lunout,*)trim(nf90_strerror(ierr)) 388 CALL ABORT_gcm("dynetat0", "", 1) 389 ENDIF 241 CALL get_var3v_t("vcov",vcov,indextime) 390 242 391 243 ! read ucov 392 corner(1)=1 393 corner(2)=1 394 corner(3)=1 395 corner(4)=indextime 396 edges(1)=iip1 397 edges(2)=jjp1 398 edges(3)=llm 399 edges(4)=1 400 ierr=nf90_inq_varid(nid,"ucov",nvarid) 401 IF (ierr .NE. nf90_noerr) THEN 402 write(lunout,*)"dynetat0: Le champ <ucov> est absent" 403 write(lunout,*)trim(nf90_strerror(ierr)) 404 CALL ABORT_gcm("dynetat0", "", 1) 405 ENDIF 406 ierr=nf90_get_var(nid,nvarid,ucov,corner,edges) 407 IF (ierr .NE. nf90_noerr) THEN 408 write(lunout,*)"dynetat0: Lecture echouee pour <ucov>" 409 write(lunout,*)trim(nf90_strerror(ierr)) 410 CALL ABORT_gcm("dynetat0", "", 1) 411 ENDIF 244 CALL get_var3u_t("ucov",ucov,indextime) 412 245 413 246 ! read teta (same corner/edges as ucov) 414 ierr=nf90_inq_varid(nid,"teta",nvarid) 415 IF (ierr .NE. nf90_noerr) THEN 416 write(lunout,*)"dynetat0: Le champ <teta> est absent" 417 write(lunout,*)trim(nf90_strerror(ierr)) 418 CALL ABORT_gcm("dynetat0", "", 1) 419 ENDIF 420 ierr=nf90_get_var(nid,nvarid,teta,corner,edges) 421 IF (ierr .NE. nf90_noerr) THEN 422 write(lunout,*)"dynetat0: Lecture echouee pour <teta>" 423 write(lunout,*)trim(nf90_strerror(ierr)) 424 CALL ABORT_gcm("dynetat0", "", 1) 425 ENDIF 247 CALL get_var3u_t("teta",teta,indextime) 426 248 427 249 ! read tracers (same corner/edges as ucov) 428 IF(nqtot.GE.1) THEN 250 corner(1)=1 251 corner(2)=1 252 corner(3)=1 253 corner(4)=indextime 254 edges(1)=iip1 255 edges(2)=jjp1 256 edges(3)=llm 257 edges(4)=1 258 IF(nqtot.GE.1) THEN 429 259 DO iq=1,nqtot 430 ierr= nf90_inq_varid( nid,tname(iq),nvarid)260 ierr= nf90_inq_varid(fID,tname(iq),vID) 431 261 IF (ierr .NE. nf90_noerr) THEN 432 write(lunout,*)"dynetat0: Le traceur <"//trim(tname(iq))// 433 & "> est absent" 434 write(lunout,*)" Il est donc initialise a zero" 262 write(lunout,*)TRIM(modname)//": Tracer <"//TRIM(var)//"> is missing" 263 write(lunout,*)" It is hence initialized to zero" 435 264 q(:,:,:,iq)=0. 265 IF (planet_type=="earth") THEN 266 !--- CRisi: for isotops, theoretical initialization using very simplified 267 ! Rayleigh distillation las. 268 IF(ok_isotopes.AND.iso_num(iq)>0) THEN 269 IF(zone_num(iq)==0) q(:,:,:,iq)=q(:,:,:,iqpere(iq))*tnat(iso_num(iq)) & 270 & *(q(:,:,:,iqpere(iq))/30.e-3)**(alpha_ideal(iso_num(iq))-1) 271 IF(zone_num(iq)==1) q(:,:,:,iq)=q(:,:,:,iqiso(iso_indnum(iq),phase_num(iq))) 272 END IF 273 ENDIF 436 274 ELSE 437 ierr=nf90_get_var( nid,nvarid,q(:,:,:,iq),corner,edges)275 ierr=nf90_get_var(fID,vID,q(:,:,:,iq),corner,edges) 438 276 IF (ierr .NE. nf90_noerr) THEN 439 write(lunout,*)"dynetat0: Lecture echouee pour " 440 &//trim(tname(iq))277 write(lunout,*)"dynetat0: Lecture echouee pour " & 278 //trim(tname(iq)) 441 279 write(lunout,*)trim(nf90_strerror(ierr)) 442 280 CALL ABORT_gcm("dynetat0", "", 1) … … 444 282 ENDIF 445 283 ENDDO 446 284 ENDIF 447 285 448 286 !read masse (same corner/edges as ucov) 449 ierr=nf90_inq_varid(nid,"masse",nvarid) 450 IF (ierr .NE. nf90_noerr) THEN 451 write(lunout,*)"dynetat0: Le champ <masse> est absent" 452 write(lunout,*)trim(nf90_strerror(ierr)) 453 CALL ABORT_gcm("dynetat0", "", 1) 454 ENDIF 455 ierr=nf90_get_var(nid,nvarid,masse,corner,edges) 456 IF (ierr .NE. nf90_noerr) THEN 457 write(lunout,*)"dynetat0: Lecture echouee pour <masse>" 458 write(lunout,*)trim(nf90_strerror(ierr)) 459 CALL ABORT_gcm("dynetat0", "", 1) 460 ENDIF 461 462 ! read ps 463 corner(1)=1 464 corner(2)=1 465 corner(3)=indextime 466 edges(1)=iip1 467 edges(2)=jjp1 468 edges(3)=1 469 ierr=nf90_inq_varid(nid,"ps",nvarid) 470 IF (ierr .NE. nf90_noerr) THEN 471 write(lunout,*)"dynetat0: Le champ <ps> est absent" 472 write(lunout,*)trim(nf90_strerror(ierr)) 473 CALL ABORT_gcm("dynetat0", "", 1) 474 ENDIF 475 ierr=nf90_get_var(nid,nvarid,ps,corner,edges) 476 IF (ierr .NE. nf90_noerr) THEN 477 write(lunout,*)"dynetat0: Lecture echouee pour <ps>" 478 write(lunout,*)trim(nf90_strerror(ierr)) 479 CALL ABORT_gcm("dynetat0", "", 1) 480 ENDIF 481 482 ierr=nf90_close(nid) 483 484 if (planet_type/="mars") then 485 day_ini=day_ini+INT(time0) ! obsolete stuff ; 0<time<1 anyways 486 time0=time0-INT(time0) 487 endif 488 489 1 FORMAT(//10x,'la valeur de im =',i4,2x,'lue sur le fichier de dem 490 *arrage est differente de la valeur parametree iim =',i4//) 491 2 FORMAT(//10x,'la valeur de jm =',i4,2x,'lue sur le fichier de dem 492 *arrage est differente de la valeur parametree jjm =',i4//) 493 3 FORMAT(//10x,'la valeur de lmax =',i4,2x,'lue sur le fichier dema 494 *rrage est differente de la valeur parametree llm =',i4//) 495 4 FORMAT(//10x,'la valeur de dtrv =',i4,2x,'lue sur le fichier dema 496 *rrage est differente de la valeur dtinteg =',i4//) 497 498 RETURN 499 END 287 CALL get_var3u_t("masse",masse,indextime) 288 289 !read ps 290 CALL get_var2_t("ps",ps,indextime) 291 292 CALL err(NF90_CLOSE(fID),"close",fichnom) 293 294 if (planet_type/="mars") then 295 day_ini=day_ini+INT(time0) ! obsolete stuff ; 0<time<1 anyways 296 time0=time0-INT(time0) 297 endif 298 299 300 CONTAINS 301 302 SUBROUTINE check_dim(n1,n2,str1,str2) 303 INTEGER, INTENT(IN) :: n1, n2 304 CHARACTER(LEN=*), INTENT(IN) :: str1, str2 305 CHARACTER(LEN=256) :: s1, s2 306 IF(n1/=n2) THEN 307 s1='value of '//TRIM(str1)//' =' 308 s2=' read in starting file differs from parametrized '//TRIM(str2)//' =' 309 WRITE(msg,'(10x,a,i4,2x,a,i4)'),s1,n1,s2,n2 310 CALL ABORT_gcm(TRIM(modname),TRIM(msg),1) 311 END IF 312 END SUBROUTINE check_dim 313 314 315 SUBROUTINE get_var1(var,v) 316 CHARACTER(LEN=*), INTENT(IN) :: var 317 REAL, INTENT(OUT) :: v(:) 318 CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) 319 CALL err(NF90_GET_VAR(fID,vID,v),"get",var) 320 END SUBROUTINE get_var1 321 322 323 SUBROUTINE get_var2(var,v) 324 CHARACTER(LEN=*), INTENT(IN) :: var 325 REAL, INTENT(OUT) :: v(:,:) 326 CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) 327 CALL err(NF90_GET_VAR(fID,vID,v),"get",var) 328 END SUBROUTINE get_var2 329 330 SUBROUTINE get_var2_t(var,v,indextime) 331 CHARACTER(LEN=*), INTENT(IN) :: var 332 REAL, INTENT(OUT) :: v(:,:) 333 INTEGER, INTENT(IN) :: indextime 334 corner(1)=1 335 corner(2)=1 336 corner(3)=indextime 337 edges(1)=iip1 338 edges(2)=jjp1 339 edges(3)=1 340 CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) 341 CALL err(NF90_GET_VAR(fID,vID,v,corner,edges),"get",var) 342 END SUBROUTINE get_var2_t 343 344 345 SUBROUTINE get_var3(var,v) ! on U grid 346 CHARACTER(LEN=*), INTENT(IN) :: var 347 REAL, INTENT(OUT) :: v(:,:,:) 348 CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) 349 CALL err(NF90_GET_VAR(fID,vID,v),"get",var) 350 END SUBROUTINE get_var3 351 352 SUBROUTINE get_var3u_t(var,v,indextime) ! on U grid 353 CHARACTER(LEN=*), INTENT(IN) :: var 354 REAL, INTENT(OUT) :: v(:,:,:) 355 INTEGER, INTENT(IN) :: indextime 356 corner(1)=1 357 corner(2)=1 358 corner(3)=1 359 corner(4)=indextime 360 edges(1)=iip1 361 edges(2)=jjp1 362 edges(3)=llm 363 edges(4)=1 364 CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) 365 CALL err(NF90_GET_VAR(fID,vID,v,corner,edges),"get",var) 366 END SUBROUTINE get_var3u_t 367 368 SUBROUTINE get_var3v_t(var,v,indextime) ! on V grid 369 CHARACTER(LEN=*), INTENT(IN) :: var 370 REAL, INTENT(OUT) :: v(:,:,:) 371 INTEGER, INTENT(IN) :: indextime 372 corner(1)=1 373 corner(2)=1 374 corner(3)=1 375 corner(4)=indextime 376 edges(1)=iip1 377 edges(2)=jjm 378 edges(3)=llm 379 edges(4)=1 380 CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var) 381 CALL err(NF90_GET_VAR(fID,vID,v,corner,edges),"get",var) 382 END SUBROUTINE get_var3v_t 383 384 SUBROUTINE err(ierr,typ,nam) 385 INTEGER, INTENT(IN) :: ierr !--- NetCDF ERROR CODE 386 CHARACTER(LEN=*), INTENT(IN) :: typ !--- TYPE OF OPERATION 387 CHARACTER(LEN=*), INTENT(IN) :: nam !--- FIELD/FILE NAME 388 IF(ierr==NF90_NoERR) RETURN 389 SELECT CASE(typ) 390 CASE('inq'); msg="Field <"//TRIM(nam)//"> is missing" 391 CASE('get'); msg="Reading failed for <"//TRIM(nam)//">" 392 CASE('open'); msg="File opening failed for <"//TRIM(nam)//">" 393 CASE('close'); msg="File closing failed for <"//TRIM(nam)//">" 394 END SELECT 395 CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr) 396 END SUBROUTINE err 397 398 END SUBROUTINE dynetat0 -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d_common/dynredem.F90 ¶
r1506 r1508 2 2 ! $Id: dynredem.F 1635 2012-07-12 11:37:16Z lguez $ 3 3 ! 4 c 5 4 ! 5 SUBROUTINE dynredem0(fichnom,iday_end,phis) 6 6 #ifdef CPP_IOIPSL 7 7 USE IOIPSL 8 8 #endif 9 USE infotrac 10 use netcdf95, only: NF95_PUT_VAR 11 use control_mod, only : planet_type 12 USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt,pa,preff, 13 . nivsig,nivsigs 14 USE comconst_mod, ONLY: daysec,dtvr,rad,omeg,g,cpp,kappa,pi 15 USE logic_mod, ONLY: fxyhypb,ysinus 16 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, 17 . taux,tauy 18 USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn,itaufin, 19 . start_time,hour_ini 20 USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 21 22 IMPLICIT NONE 23 c======================================================================= 24 c Ecriture du fichier de redemarrage sous format NetCDF (initialisation) 25 c======================================================================= 26 c Declarations: 27 c ------------- 28 #include "dimensions.h" 29 #include "paramet.h" 30 #include "comgeom2.h" 31 #include "netcdf.inc" 32 #include "iniprint.h" 33 34 c Arguments: 35 c ---------- 36 INTEGER iday_end 37 REAL phis(iip1, jjp1) 38 CHARACTER*(*) fichnom 39 40 c Local: 41 c ------ 42 INTEGER iq,l 43 INTEGER length 44 PARAMETER (length = 100) 45 REAL tab_cntrl(length) ! tableau des parametres du run 46 INTEGER ierr 47 character*20 modname 48 character*80 abort_message 49 50 c Variables locales pour NetCDF: 51 c 52 INTEGER dims2(2), dims3(3), dims4(4) 53 INTEGER idim_index 54 INTEGER idim_rlonu, idim_rlonv, idim_rlatu, idim_rlatv 55 INTEGER idim_s, idim_sig 56 INTEGER idim_tim 57 INTEGER nid,nvarid 58 59 REAL zan0,zjulian,hours 60 INTEGER yyears0,jjour0, mmois0 61 character*30 unites 62 63 character(len=12) :: start_file_type="earth" ! default start file type 64 INTEGER idecal 65 66 c----------------------------------------------------------------------- 67 modname='dynredem0' 9 USE infotrac, ONLY: nqtot, tname, ttext 10 USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL, & 11 NF90_CLOSE, NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER 12 USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil 13 use netcdf95, only: NF95_PUT_VAR 14 use control_mod, only : planet_type 15 USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt,pa,preff, & 16 nivsig,nivsigs 17 USE comconst_mod, ONLY: daysec,dtvr,rad,omeg,g,cpp,kappa,pi 18 USE logic_mod, ONLY: fxyhypb,ysinus 19 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, & 20 taux,tauy 21 USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn,itaufin, & 22 start_time,hour_ini 23 USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 24 25 IMPLICIT NONE 26 !======================================================================= 27 ! Writting the NetCDF restart file (initialisation) 28 !======================================================================= 29 ! Declarations: 30 ! ------------- 31 include "dimensions.h" 32 include "paramet.h" 33 include "comgeom2.h" 34 include "netcdf.inc" 35 include "iniprint.h" 36 37 !=============================================================================== 38 ! Arguments: 39 CHARACTER(LEN=*), INTENT(IN) :: fichnom !--- FILE NAME 40 INTEGER, INTENT(IN) :: iday_end !--- 41 REAL, INTENT(IN) :: phis(iip1, jjp1) !--- GROUND GEOPOTENTIAL 42 !=============================================================================== 43 ! Local variables: 44 INTEGER :: iq,l 45 INTEGER, PARAMETER :: length=100 46 REAL :: tab_cntrl(length) ! run parameters 47 INTEGER :: ierr 48 CHARACTER(LEN=80) :: abort_message 49 50 ! For NetCDF: 51 CHARACTER(LEN=30) :: unites 52 INTEGER :: indexID 53 INTEGER :: rlonuID, rlonvID, rlatuID, rlatvID 54 INTEGER :: sID, sigID, nID, vID, timID 55 INTEGER :: yyears0, jjour0, mmois0 56 REAL :: zan0, zjulian, hours 57 58 CHARACTER(len=12) :: start_file_type="earth" ! default start file type 59 INTEGER :: idecal 60 61 !=============================================================================== 62 ! fill dynredem_mod module variables 63 modname='dynredem0'; fil=fichnom 68 64 69 65 #ifdef CPP_IOIPSL 70 71 66 call ymds2ju(annee_ref, 1, iday_end, 0.0, zjulian) 67 call ju2ymds(zjulian, yyears0, mmois0, jjour0, hours) 72 68 #else 73 69 ! set yyears0, mmois0, jjour0 to 0,1,1 (hours is not used) 74 75 76 70 yyears0=0 71 mmois0=1 72 jjour0=1 77 73 #endif 78 74 79 !!! AS: idecal is a hack to be able to read planeto starts... 80 !!! .... while keeping everything OK for LMDZ EARTH 81 if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then 82 write(lunout,*) trim(modname),' : Planeto-like start file' 83 start_file_type="planeto" 84 idecal = 4 85 else 86 write(lunout,*) trim(modname),' : Earth-like start file' 87 idecal = 5 88 endif 89 90 DO l=1,length 91 tab_cntrl(l) = 0. 92 ENDDO 93 tab_cntrl(1) = REAL(iim) 94 tab_cntrl(2) = REAL(jjm) 95 tab_cntrl(3) = REAL(llm) 96 if (start_file_type.eq."earth") then 97 tab_cntrl(4)=REAL(day_ref) 98 else 99 !tab_cntrl(4)=REAL(day_end) 100 tab_cntrl(4)=REAL(iday_end) 101 endif 102 tab_cntrl(5) = REAL(annee_ref) 103 tab_cntrl(idecal+1) = rad 104 tab_cntrl(idecal+2) = omeg 105 tab_cntrl(idecal+3) = g 106 tab_cntrl(idecal+4) = cpp 107 tab_cntrl(idecal+5) = kappa 108 tab_cntrl(idecal+6) = daysec 109 tab_cntrl(idecal+7) = dtvr 110 tab_cntrl(idecal+8) = etot0 111 tab_cntrl(idecal+9) = ptot0 112 tab_cntrl(idecal+10) = ztot0 113 tab_cntrl(idecal+11) = stot0 114 tab_cntrl(idecal+12) = ang0 115 tab_cntrl(idecal+13) = pa 116 tab_cntrl(idecal+14) = preff 117 c 118 c ..... parametres pour le zoom ...... 119 120 tab_cntrl(idecal+15) = clon 121 tab_cntrl(idecal+16) = clat 122 tab_cntrl(idecal+17) = grossismx 123 tab_cntrl(idecal+18) = grossismy 124 c 125 IF ( fxyhypb ) THEN 126 tab_cntrl(idecal+19) = 1. 127 tab_cntrl(idecal+20) = dzoomx 128 tab_cntrl(idecal+21) = dzoomy 129 tab_cntrl(idecal+22) = 0. 130 tab_cntrl(idecal+23) = taux 131 tab_cntrl(idecal+24) = tauy 132 ELSE 133 tab_cntrl(idecal+19) = 0. 134 tab_cntrl(idecal+20) = dzoomx 135 tab_cntrl(idecal+21) = dzoomy 136 tab_cntrl(idecal+22) = 0. 137 tab_cntrl(idecal+23) = 0. 138 tab_cntrl(idecal+24) = 0. 139 IF( ysinus ) tab_cntrl(idecal+22) = 1. 140 ENDIF 141 142 if (start_file_type.eq."earth") then 143 tab_cntrl(idecal+25) = REAL(iday_end) 144 tab_cntrl(idecal+26) = REAL(itau_dyn + itaufin) 145 c start_time: start_time of simulation (not necessarily 0.) 146 tab_cntrl(idecal+27) = start_time 147 endif 148 149 if (planet_type=="mars") then ! For Mars only 150 tab_cntrl(29)=hour_ini 151 endif 152 c 153 c ......................................................... 154 c 155 c Creation du fichier: 156 c 157 ierr = NF_CREATE(fichnom, NF_CLOBBER, nid) 158 IF (ierr.NE.NF_NOERR) THEN 159 write(lunout,*)"dynredem0: Pb d ouverture du fichier " 160 & //trim(fichnom) 161 write(lunout,*)' ierr = ', ierr 162 CALL ABORT_GCM("DYNREDEM0", "", 1) 163 ENDIF 164 c 165 c Preciser quelques attributs globaux: 166 c 167 ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 27, 168 . "Fichier demmarage dynamique") 169 c 170 c Definir les dimensions du fichiers: 171 c 172 if (start_file_type.eq."earth") then 173 ierr = NF_DEF_DIM (nid, "index", length, idim_index) 174 ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu) 175 ierr = NF_DEF_DIM (nid, "rlatu", jjp1, idim_rlatu) 176 ierr = NF_DEF_DIM (nid, "rlonv", iip1, idim_rlonv) 177 ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv) 178 ierr = NF_DEF_DIM (nid, "sigs", llm, idim_s) 179 ierr = NF_DEF_DIM (nid, "sig", llmp1, idim_sig) 180 ierr = NF_DEF_DIM (nid, "temps", NF_UNLIMITED, idim_tim) 181 else 182 ierr = NF_DEF_DIM (nid, "index", length, idim_index) 183 ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu) 184 ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu) 185 ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv) 186 ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv) 187 ierr = NF_DEF_DIM (nid, "altitude", llm, idim_s) 188 ierr = NF_DEF_DIM (nid, "interlayer", llmp1, idim_sig) 189 ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_tim) 190 endif 191 c 192 ierr = NF_ENDDEF(nid) ! sortir du mode de definition 193 c 194 c Definir et enregistrer certains champs invariants: 195 c 196 ierr = NF_REDEF (nid) 197 cIM 220306 BEG 198 #ifdef NC_DOUBLE 199 ierr = NF_DEF_VAR (nid,"controle",NF_DOUBLE,1,idim_index,nvarid) 200 #else 201 ierr = NF_DEF_VAR (nid,"controle",NF_FLOAT,1,idim_index,nvarid) 202 #endif 203 cIM 220306 END 204 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, 205 . "Parametres de controle") 206 ierr = NF_ENDDEF(nid) 207 call NF95_PUT_VAR(nid,nvarid,tab_cntrl) 208 c 209 ierr = NF_REDEF (nid) 210 cIM 220306 BEG 211 #ifdef NC_DOUBLE 212 ierr = NF_DEF_VAR (nid,"rlonu",NF_DOUBLE,1,idim_rlonu,nvarid) 213 #else 214 ierr = NF_DEF_VAR (nid,"rlonu",NF_FLOAT,1,idim_rlonu,nvarid) 215 #endif 216 cIM 220306 END 217 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23, 218 . "Longitudes des points U") 219 ierr = NF_ENDDEF(nid) 220 call NF95_PUT_VAR(nid,nvarid,rlonu) 221 c 222 ierr = NF_REDEF (nid) 223 cIM 220306 BEG 224 #ifdef NC_DOUBLE 225 ierr = NF_DEF_VAR (nid,"rlatu",NF_DOUBLE,1,idim_rlatu,nvarid) 226 #else 227 ierr = NF_DEF_VAR (nid,"rlatu",NF_FLOAT,1,idim_rlatu,nvarid) 228 #endif 229 cIM 220306 END 230 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, 231 . "Latitudes des points U") 232 ierr = NF_ENDDEF(nid) 233 call NF95_PUT_VAR (nid,nvarid,rlatu) 234 c 235 ierr = NF_REDEF (nid) 236 cIM 220306 BEG 237 #ifdef NC_DOUBLE 238 ierr = NF_DEF_VAR (nid,"rlonv",NF_DOUBLE,1,idim_rlonv,nvarid) 239 #else 240 ierr = NF_DEF_VAR (nid,"rlonv",NF_FLOAT,1,idim_rlonv,nvarid) 241 #endif 242 cIM 220306 END 243 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23, 244 . "Longitudes des points V") 245 ierr = NF_ENDDEF(nid) 246 call NF95_PUT_VAR(nid,nvarid,rlonv) 247 c 248 ierr = NF_REDEF (nid) 249 cIM 220306 BEG 250 #ifdef NC_DOUBLE 251 ierr = NF_DEF_VAR (nid,"rlatv",NF_DOUBLE,1,idim_rlatv,nvarid) 252 #else 253 ierr = NF_DEF_VAR (nid,"rlatv",NF_FLOAT,1,idim_rlatv,nvarid) 254 #endif 255 cIM 220306 END 256 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, 257 . "Latitudes des points V") 258 ierr = NF_ENDDEF(nid) 259 call NF95_PUT_VAR(nid,nvarid,rlatv) 260 c 261 if (start_file_type.eq."earth") then 262 ierr = NF_REDEF (nid) 263 cIM 220306 BEG 264 #ifdef NC_DOUBLE 265 ierr = NF_DEF_VAR (nid,"nivsigs",NF_DOUBLE,1,idim_s,nvarid) 266 #else 267 ierr = NF_DEF_VAR (nid,"nivsigs",NF_FLOAT,1,idim_s,nvarid) 268 #endif 269 cIM 220306 END 270 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 28, 271 . "Numero naturel des couches s") 272 ierr = NF_ENDDEF(nid) 273 call NF95_PUT_VAR(nid,nvarid,nivsigs) 274 c 275 ierr = NF_REDEF (nid) 276 cIM 220306 BEG 277 #ifdef NC_DOUBLE 278 ierr = NF_DEF_VAR (nid,"nivsig",NF_DOUBLE,1,idim_sig,nvarid) 279 #else 280 ierr = NF_DEF_VAR (nid,"nivsig",NF_FLOAT,1,idim_sig,nvarid) 281 #endif 282 cIM 220306 END 283 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 32, 284 . "Numero naturel des couches sigma") 285 ierr = NF_ENDDEF(nid) 286 call NF95_PUT_VAR(nid,nvarid,nivsig) 287 endif ! of if (start_file_type.eq."earth") 288 c 289 ierr = NF_REDEF (nid) 290 cIM 220306 BEG 291 #ifdef NC_DOUBLE 292 ierr = NF_DEF_VAR (nid,"ap",NF_DOUBLE,1,idim_sig,nvarid) 293 #else 294 ierr = NF_DEF_VAR (nid,"ap",NF_FLOAT,1,idim_sig,nvarid) 295 #endif 296 cIM 220306 END 297 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26, 298 . "Coefficient A pour hybride") 299 ierr = NF_ENDDEF(nid) 300 call NF95_PUT_VAR(nid,nvarid,ap) 301 c 302 ierr = NF_REDEF (nid) 303 cIM 220306 BEG 304 #ifdef NC_DOUBLE 305 ierr = NF_DEF_VAR (nid,"bp",NF_DOUBLE,1,idim_sig,nvarid) 306 #else 307 ierr = NF_DEF_VAR (nid,"bp",NF_FLOAT,1,idim_sig,nvarid) 308 #endif 309 cIM 220306 END 310 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26, 311 . "Coefficient B pour hybride") 312 ierr = NF_ENDDEF(nid) 313 call NF95_PUT_VAR(nid,nvarid,bp) 314 c 315 if (start_file_type.ne."earth") then 316 ierr = NF_REDEF (nid) 317 cIM 220306 BEG 318 #ifdef NC_DOUBLE 319 ierr = NF_DEF_VAR (nid,"aps",NF_DOUBLE,1,idim_s,nvarid) 320 #else 321 ierr = NF_DEF_VAR (nid,"aps",NF_FLOAT,1,idim_s,nvarid) 322 #endif 323 cIM 220306 END 324 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 37, 325 . "Coef AS: hybrid pressure at midlayers") 326 ierr = NF_ENDDEF(nid) 327 call NF95_PUT_VAR(nid,nvarid,aps) 328 c 329 ierr = NF_REDEF (nid) 330 cIM 220306 BEG 331 #ifdef NC_DOUBLE 332 ierr = NF_DEF_VAR (nid,"bps",NF_DOUBLE,1,idim_s,nvarid) 333 #else 334 ierr = NF_DEF_VAR (nid,"bps",NF_FLOAT,1,idim_s,nvarid) 335 #endif 336 cIM 220306 END 337 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 34, 338 . "Coef BS: hybrid sigma at midlayers") 339 ierr = NF_ENDDEF(nid) 340 call NF95_PUT_VAR(nid,nvarid,bps) 341 endif ! of if (start_file_type.ne."earth") 342 c 343 ierr = NF_REDEF (nid) 344 cIM 220306 BEG 345 #ifdef NC_DOUBLE 346 ierr = NF_DEF_VAR (nid,"presnivs",NF_DOUBLE,1,idim_s,nvarid) 347 #else 348 ierr = NF_DEF_VAR (nid,"presnivs",NF_FLOAT,1,idim_s,nvarid) 349 #endif 350 cIM 220306 END 351 ierr = NF_ENDDEF(nid) 352 call NF95_PUT_VAR(nid,nvarid,presnivs) 353 c 354 if (start_file_type.ne."earth") then 75 !!! AS: idecal is a hack to be able to read planeto starts... 76 !!! .... while keeping everything OK for LMDZ EARTH 77 if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then 78 write(lunout,*) trim(modname),' : Planeto-like start file' 79 start_file_type="planeto" 80 idecal = 4 81 else 82 write(lunout,*) trim(modname),' : Earth-like start file' 83 idecal = 5 84 endif 85 86 tab_cntrl(:) = 0. 87 tab_cntrl(1) = REAL(iim) 88 tab_cntrl(2) = REAL(jjm) 89 tab_cntrl(3) = REAL(llm) 90 if (start_file_type.eq."earth") then 91 tab_cntrl(4)=REAL(day_ref) 92 else 93 !tab_cntrl(4)=REAL(day_end) 94 tab_cntrl(4)=REAL(iday_end) 95 endif 96 tab_cntrl(5) = REAL(annee_ref) 97 tab_cntrl(idecal+1) = rad 98 tab_cntrl(idecal+2) = omeg 99 tab_cntrl(idecal+3) = g 100 tab_cntrl(idecal+4) = cpp 101 tab_cntrl(idecal+5) = kappa 102 tab_cntrl(idecal+6) = daysec 103 tab_cntrl(idecal+7) = dtvr 104 tab_cntrl(idecal+8) = etot0 105 tab_cntrl(idecal+9) = ptot0 106 tab_cntrl(idecal+10) = ztot0 107 tab_cntrl(idecal+11) = stot0 108 tab_cntrl(idecal+12) = ang0 109 tab_cntrl(idecal+13) = pa 110 tab_cntrl(idecal+14) = preff 111 112 ! ..... parameters for the zoom ...... 113 tab_cntrl(idecal+15) = clon 114 tab_cntrl(idecal+16) = clat 115 tab_cntrl(idecal+17) = grossismx 116 tab_cntrl(idecal+18) = grossismy 117 ! 118 IF ( fxyhypb ) THEN 119 tab_cntrl(idecal+19) = 1. 120 tab_cntrl(idecal+20) = dzoomx 121 tab_cntrl(idecal+21) = dzoomy 122 tab_cntrl(idecal+22) = 0. 123 tab_cntrl(idecal+23) = taux 124 tab_cntrl(idecal+24) = tauy 125 ELSE 126 tab_cntrl(idecal+19) = 0. 127 tab_cntrl(idecal+20) = dzoomx 128 tab_cntrl(idecal+21) = dzoomy 129 tab_cntrl(idecal+22) = 0. 130 tab_cntrl(idecal+23) = 0. 131 tab_cntrl(idecal+24) = 0. 132 IF( ysinus ) tab_cntrl(idecal+22) = 1. 133 ENDIF 134 135 if (start_file_type.eq."earth") then 136 tab_cntrl(idecal+25) = REAL(iday_end) 137 tab_cntrl(idecal+26) = REAL(itau_dyn + itaufin) 138 ! start_time: start_time of simulation (not necessarily 0.) 139 tab_cntrl(idecal+27) = start_time 140 endif 141 142 if (planet_type=="mars") then ! For Mars only 143 tab_cntrl(29)=hour_ini 144 endif 145 146 !--- File creation 147 CALL err(NF90_CREATE(fichnom,NF90_CLOBBER,nid)) 148 149 !--- Some global attributes 150 CALL err(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier demarrage dynamique")) 151 152 !--- Dimensions 153 if (start_file_type.eq."earth") then 154 CALL err(NF90_DEF_DIM(nid,"index", length, indexID)) 155 CALL err(NF90_DEF_DIM(nid,"rlonu", iip1, rlonuID)) 156 CALL err(NF90_DEF_DIM(nid,"rlatu", jjp1, rlatuID)) 157 CALL err(NF90_DEF_DIM(nid,"rlonv", iip1, rlonvID)) 158 CALL err(NF90_DEF_DIM(nid,"rlatv", jjm, rlatvID)) 159 CALL err(NF90_DEF_DIM(nid,"sigs", llm, sID)) 160 CALL err(NF90_DEF_DIM(nid,"sig", llmp1, sigID)) 161 CALL err(NF90_DEF_DIM(nid,"temps", NF90_UNLIMITED, timID)) 162 else 163 CALL err(NF90_DEF_DIM(nid,"index", length, indexID)) 164 CALL err(NF90_DEF_DIM(nid,"rlonu", iip1, rlonuID)) 165 CALL err(NF90_DEF_DIM(nid,"latitude", jjp1, rlatuID)) 166 CALL err(NF90_DEF_DIM(nid,"longitude", iip1, rlonvID)) 167 CALL err(NF90_DEF_DIM(nid,"rlatv", jjm, rlatvID)) 168 CALL err(NF90_DEF_DIM(nid,"altitude", llm, sID)) 169 CALL err(NF90_DEF_DIM(nid,"interlayer", llmp1, sigID)) 170 CALL err(NF90_DEF_DIM(nid,"Time", NF90_UNLIMITED, timID)) 171 endif 172 173 !--- Define and save invariant fields 174 CALL put_var1(nid,"controle","Parametres de controle" ,[indexID],tab_cntrl) 175 CALL put_var1(nid,"rlonu" ,"Longitudes des points U",[rlonuID],rlonu) 176 CALL put_var1(nid,"rlatu" ,"Latitudes des points U" ,[rlatuID],rlatu) 177 CALL put_var1(nid,"rlonv" ,"Longitudes des points V",[rlonvID],rlonv) 178 CALL put_var1(nid,"rlatv" ,"Latitudes des points V" ,[rlatvID],rlatv) 179 if (start_file_type.eq."earth") then 180 CALL put_var1(nid,"nivsigs" ,"Numero naturel des couches s" ,[sID] ,nivsigs) 181 CALL put_var1(nid,"nivsig" ,"Numero naturel des couches sigma",[sigID],nivsig) 182 endif ! of if (start_file_type.eq."earth") 183 CALL put_var1(nid,"ap" ,"Coefficient A pour hybride" ,[sigID],ap) 184 CALL put_var1(nid,"bp" ,"Coefficient B pour hybride" ,[sigID],bp) 185 if (start_file_type.ne."earth") then 186 CALL put_var1(nid,"aps","Coef AS: hybrid pressure at midlayers",[sID],aps) 187 CALL put_var1(nid,"bps","Coef BS: hybrid sigma at midlayers",[sID],bps) 188 endif ! of if (start_file_type.eq."earth") 189 CALL put_var1(nid,"presnivs","" ,[sID] ,presnivs) 190 if (start_file_type.ne."earth") then 355 191 ierr = NF_REDEF (nid) 356 192 #ifdef NC_DOUBLE 357 ierr = NF_DEF_VAR(nid,"latitude",NF_DOUBLE,1, idim_rlatu,nvarid)193 ierr = NF_DEF_VAR(nid,"latitude",NF_DOUBLE,1,rlatuID,vID) 358 194 #else 359 ierr = NF_DEF_VAR(nid,"latitude",NF_FLOAT,1, idim_rlatu,nvarid)195 ierr = NF_DEF_VAR(nid,"latitude",NF_FLOAT,1,rlatuID,vID) 360 196 #endif 361 ierr =NF_PUT_ATT_TEXT(nid, nvarid,'units',13,"degrees_north")362 ierr = NF_PUT_ATT_TEXT (nid, nvarid,"long_name", 14,363 ."North latitude")197 ierr =NF_PUT_ATT_TEXT(nid,vID,'units',13,"degrees_north") 198 ierr = NF_PUT_ATT_TEXT (nid,vID,"long_name", 14, & 199 "North latitude") 364 200 ierr = NF_ENDDEF(nid) 365 call NF95_PUT_VAR(nid, nvarid,rlatu*180/pi)366 c 201 call NF95_PUT_VAR(nid,vID,rlatu*180/pi) 202 ! 367 203 ierr = NF_REDEF (nid) 368 204 #ifdef NC_DOUBLE 369 ierr=NF_DEF_VAR(nid,"longitude",NF_DOUBLE,1, idim_rlonv,nvarid)205 ierr=NF_DEF_VAR(nid,"longitude",NF_DOUBLE,1,rlonvID,vID) 370 206 #else 371 ierr=NF_DEF_VAR(nid,"longitude",NF_FLOAT,1, idim_rlonv,nvarid)207 ierr=NF_DEF_VAR(nid,"longitude",NF_FLOAT,1,rlonvID,vID) 372 208 #endif 373 ierr = NF_PUT_ATT_TEXT (nid, nvarid,"long_name", 14,374 ."East longitude")375 ierr = NF_PUT_ATT_TEXT(nid, nvarid,'units',12,"degrees_east")209 ierr = NF_PUT_ATT_TEXT (nid,vID,"long_name", 14, & 210 "East longitude") 211 ierr = NF_PUT_ATT_TEXT(nid,vID,'units',12,"degrees_east") 376 212 ierr = NF_ENDDEF(nid) 377 call NF95_PUT_VAR(nid, nvarid,rlonv*180/pi)378 c 213 call NF95_PUT_VAR(nid,vID,rlonv*180/pi) 214 ! 379 215 ierr = NF_REDEF (nid) 380 216 #ifdef NC_DOUBLE 381 ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, 382 . idim_s,nvarid)217 ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, & 218 sID,vID) 383 219 #else 384 ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, 385 . idim_s,nvarid)220 ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, & 221 sID,vID) 386 222 #endif 387 ierr = NF_PUT_ATT_TEXT(nid, nvarid,"long_name",10,"pseudo-alt")388 ierr = NF_PUT_ATT_TEXT (nid, nvarid,'units',2,"km")389 ierr = NF_PUT_ATT_TEXT (nid, nvarid,'positive',2,"up")223 ierr = NF_PUT_ATT_TEXT(nid,vID,"long_name",10,"pseudo-alt") 224 ierr = NF_PUT_ATT_TEXT (nid,vID,'units',2,"km") 225 ierr = NF_PUT_ATT_TEXT (nid,vID,'positive',2,"up") 390 226 ierr = NF_ENDDEF(nid) 391 call NF95_PUT_VAR(nid,nvarid,pseudoalt) 392 endif ! of if (start_file_type.ne."earth") 393 c 394 c Coefficients de passage cov. <-> contra. <--> naturel 395 c 396 ierr = NF_REDEF (nid) 397 dims2(1) = idim_rlonu 398 dims2(2) = idim_rlatu 399 cIM 220306 BEG 400 #ifdef NC_DOUBLE 401 ierr = NF_DEF_VAR (nid,"cu",NF_DOUBLE,2,dims2,nvarid) 402 #else 403 ierr = NF_DEF_VAR (nid,"cu",NF_FLOAT,2,dims2,nvarid) 404 #endif 405 cIM 220306 END 406 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29, 407 . "Coefficient de passage pour U") 408 ierr = NF_ENDDEF(nid) 409 call NF95_PUT_VAR(nid,nvarid,cu) 410 c 411 ierr = NF_REDEF (nid) 412 dims2(1) = idim_rlonv 413 dims2(2) = idim_rlatv 414 cIM 220306 BEG 415 #ifdef NC_DOUBLE 416 ierr = NF_DEF_VAR (nid,"cv",NF_DOUBLE,2,dims2,nvarid) 417 #else 418 ierr = NF_DEF_VAR (nid,"cv",NF_FLOAT,2,dims2,nvarid) 419 #endif 420 cIM 220306 END 421 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29, 422 . "Coefficient de passage pour V") 423 ierr = NF_ENDDEF(nid) 424 call NF95_PUT_VAR(nid,nvarid,cv) 425 c 426 c Aire de chaque maille: 427 c 428 ierr = NF_REDEF (nid) 429 dims2(1) = idim_rlonv 430 dims2(2) = idim_rlatu 431 cIM 220306 BEG 432 #ifdef NC_DOUBLE 433 ierr = NF_DEF_VAR (nid,"aire",NF_DOUBLE,2,dims2,nvarid) 434 #else 435 ierr = NF_DEF_VAR (nid,"aire",NF_FLOAT,2,dims2,nvarid) 436 #endif 437 cIM 220306 END 438 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, 439 . "Aires de chaque maille") 440 ierr = NF_ENDDEF(nid) 441 call NF95_PUT_VAR(nid,nvarid,aire) 442 c 443 c Geopentiel au sol: 444 c 445 ierr = NF_REDEF (nid) 446 dims2(1) = idim_rlonv 447 dims2(2) = idim_rlatu 448 cIM 220306 BEG 449 #ifdef NC_DOUBLE 450 ierr = NF_DEF_VAR (nid,"phisinit",NF_DOUBLE,2,dims2,nvarid) 451 #else 452 ierr = NF_DEF_VAR (nid,"phisinit",NF_FLOAT,2,dims2,nvarid) 453 #endif 454 cIM 220306 END 455 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19, 456 . "Geopotentiel au sol") 457 ierr = NF_ENDDEF(nid) 458 call NF95_PUT_VAR(nid,nvarid,phis) 459 c 460 c Definir les variables pour pouvoir les enregistrer plus tard: 461 c 462 ierr = NF_REDEF (nid) ! entrer dans le mode de definition 463 c 464 if (start_file_type.eq."earth") then 465 cIM 220306 BEG 466 #ifdef NC_DOUBLE 467 ierr = NF_DEF_VAR (nid,"temps",NF_DOUBLE,1,idim_tim,nvarid) 468 #else 469 ierr = NF_DEF_VAR (nid,"temps",NF_FLOAT,1,idim_tim,nvarid) 470 #endif 471 cIM 220306 END 472 else ! start_file_type=="planeto" 473 #ifdef NC_DOUBLE 474 ierr = NF_DEF_VAR (nid,"Time",NF_DOUBLE,1,idim_tim,nvarid) 475 #else 476 ierr = NF_DEF_VAR (nid,"Time",NF_FLOAT,1,idim_tim,nvarid) 477 #endif 478 endif ! of if (start_file_type.eq."earth") 479 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19, 480 . "Temps de simulation") 481 write(unites,200)yyears0,mmois0,jjour0 482 200 format('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00') 483 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "units", 30, 484 . unites) 485 486 c 487 dims4(1) = idim_rlonu 488 dims4(2) = idim_rlatu 489 dims4(3) = idim_s 490 dims4(4) = idim_tim 491 cIM 220306 BEG 492 #ifdef NC_DOUBLE 493 ierr = NF_DEF_VAR (nid,"ucov",NF_DOUBLE,4,dims4,nvarid) 494 #else 495 ierr = NF_DEF_VAR (nid,"ucov",NF_FLOAT,4,dims4,nvarid) 496 #endif 497 cIM 220306 END 498 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9, 499 . "Vitesse U") 500 c 501 dims4(1) = idim_rlonv 502 dims4(2) = idim_rlatv 503 dims4(3) = idim_s 504 dims4(4) = idim_tim 505 cIM 220306 BEG 506 #ifdef NC_DOUBLE 507 ierr = NF_DEF_VAR (nid,"vcov",NF_DOUBLE,4,dims4,nvarid) 508 #else 509 ierr = NF_DEF_VAR (nid,"vcov",NF_FLOAT,4,dims4,nvarid) 510 #endif 511 cIM 220306 END 512 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9, 513 . "Vitesse V") 514 c 515 dims4(1) = idim_rlonv 516 dims4(2) = idim_rlatu 517 dims4(3) = idim_s 518 dims4(4) = idim_tim 519 cIM 220306 BEG 520 #ifdef NC_DOUBLE 521 ierr = NF_DEF_VAR (nid,"teta",NF_DOUBLE,4,dims4,nvarid) 522 #else 523 ierr = NF_DEF_VAR (nid,"teta",NF_FLOAT,4,dims4,nvarid) 524 #endif 525 cIM 220306 END 526 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 11, 527 . "Temperature") 528 c 529 dims4(1) = idim_rlonv 530 dims4(2) = idim_rlatu 531 dims4(3) = idim_s 532 dims4(4) = idim_tim 533 IF(nqtot.GE.1) THEN 534 DO iq=1,nqtot 535 cIM 220306 BEG 536 #ifdef NC_DOUBLE 537 ierr = NF_DEF_VAR (nid,tname(iq),NF_DOUBLE,4,dims4,nvarid) 538 #else 539 ierr = NF_DEF_VAR (nid,tname(iq),NF_FLOAT,4,dims4,nvarid) 540 #endif 541 cIM 220306 END 542 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,ttext(iq)) 543 ENDDO 544 ENDIF 545 c 546 dims4(1) = idim_rlonv 547 dims4(2) = idim_rlatu 548 dims4(3) = idim_s 549 dims4(4) = idim_tim 550 cIM 220306 BEG 551 #ifdef NC_DOUBLE 552 ierr = NF_DEF_VAR (nid,"masse",NF_DOUBLE,4,dims4,nvarid) 553 #else 554 ierr = NF_DEF_VAR (nid,"masse",NF_FLOAT,4,dims4,nvarid) 555 #endif 556 cIM 220306 END 557 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12, 558 . "C est quoi ?") 559 c 560 dims3(1) = idim_rlonv 561 dims3(2) = idim_rlatu 562 dims3(3) = idim_tim 563 cIM 220306 BEG 564 #ifdef NC_DOUBLE 565 ierr = NF_DEF_VAR (nid,"ps",NF_DOUBLE,3,dims3,nvarid) 566 #else 567 ierr = NF_DEF_VAR (nid,"ps",NF_FLOAT,3,dims3,nvarid) 568 #endif 569 cIM 220306 END 570 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 15, 571 . "Pression au sol") 572 c 573 ierr = NF_ENDDEF(nid) ! sortir du mode de definition 574 ierr = NF_CLOSE(nid) ! fermer le fichier 575 576 write(lunout,*)'dynredem0: iim,jjm,llm,iday_end', 577 & iim,jjm,llm,iday_end 578 write(lunout,*)'dynredem0: rad,omeg,g,cpp,kappa', 579 & rad,omeg,g,cpp,kappa 580 581 RETURN 582 END 227 call NF95_PUT_VAR(nid,vID,pseudoalt) 228 CALL err(NF_REDEF(nid)) 229 endif ! of if (start_file_type.ne."earth") 230 231 ! covariant <-> contravariant <-> natural conversion coefficients 232 CALL put_var2(nid,"cu","Coefficient de passage pour U",[rlonuID,rlatuID],cu) 233 CALL put_var2(nid,"cv","Coefficient de passage pour V",[rlonvID,rlatvID],cv) 234 CALL put_var2(nid,"aire","Aires de chaque maille" ,[rlonvID,rlatuID],aire) 235 CALL put_var2(nid,"phisinit","Geopotentiel au sol" ,[rlonvID,rlatuID],phis) 236 237 238 ! Define variables that will be stored later: 239 WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')"),& 240 yyears0,mmois0,jjour0 241 IF (planet_type.eq."earth") THEN 242 CALL cre_var(nid,"temps","Temps de simulation",[timID],unites) 243 ELSE 244 CALL cre_var(nid,"Time","Temps de simulation",[timID],unites) 245 ENDIF 246 247 CALL cre_var(nid,"ucov" ,"Vitesse U" ,[rlonuID,rlatuID,sID,timID]) 248 CALL cre_var(nid,"vcov" ,"Vitesse V" ,[rlonvID,rlatvID,sID,timID]) 249 CALL cre_var(nid,"teta" ,"Temperature",[rlonvID,rlatuID,sID,timID]) 250 251 IF(nqtot.GE.1) THEN 252 DO iq=1,nqtot 253 CALL cre_var(nid,tname(iq),ttext(iq),[rlonvID,rlatuID,sID,timID]) 254 END DO 255 ENDIF 256 257 CALL cre_var(nid,"masse","Masse d air" ,[rlonvID,rlatuID,sID,timID]) 258 CALL cre_var(nid,"ps" ,"Pression au sol",[rlonvID,rlatuID ,timID]) 259 260 CALL err(NF90_CLOSE (nid)) ! close file 261 262 WRITE(lunout,*)TRIM(modname)//': iim,jjm,llm,iday_end',iim,jjm,llm,iday_end 263 WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa 264 265 END SUBROUTINE dynredem0 583 266 584 267 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 585 268 586 SUBROUTINE dynredem1(fichnom,time, 587 . vcov,ucov,teta,q,masse,ps) 588 USE infotrac 589 USE control_mod, only : planet_type 590 use netcdf, only: NF90_get_VAR 591 use netcdf95, only: NF95_PUT_VAR 592 USE temps_mod, ONLY: itaufin,itau_dyn 269 SUBROUTINE dynredem1(fichnom,time,vcov,ucov,teta,q,masse,ps) 270 ! 271 !------------------------------------------------------------------------------- 272 ! Purpose: Write the NetCDF restart file (append). 273 !------------------------------------------------------------------------------- 274 USE infotrac, ONLY: nqtot, tname, type_trac 275 USE control_mod, only : planet_type 276 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_GET_VAR, NF90_INQ_VARID, & 277 NF90_CLOSE, NF90_WRITE, NF90_PUT_VAR, NF90_NoErr 278 use netcdf95, only: NF95_PUT_VAR 279 USE temps_mod, ONLY: itaufin,itau_dyn 280 USE dynredem_mod, ONLY: dynredem_write_u, dynredem_write_v, dynredem_read_u, & 281 err, modname, fil, msg 593 282 594 IMPLICIT NONE 595 c================================================================= 596 c Ecriture du fichier de redemarrage sous format NetCDF 597 c================================================================= 598 #include "dimensions.h" 599 #include "paramet.h" 600 #include "netcdf.inc" 601 #include "comgeom.h" 602 #include "iniprint.h" 603 604 605 INTEGER l 606 REAL vcov(iip1,jjm,llm),ucov(iip1, jjp1,llm) 607 REAL teta(iip1, jjp1,llm) 608 REAL ps(iip1, jjp1),masse(iip1, jjp1,llm) 609 REAL q(iip1, jjp1, llm, nqtot) 610 CHARACTER*(*) fichnom 611 612 REAL time 613 INTEGER nid, nvarid, nid_trac, nvarid_trac 614 REAL trac_tmp(ip1jmp1,llm) 615 INTEGER ierr, ierr_file 616 INTEGER iq 617 INTEGER length 618 PARAMETER (length = 100) 619 REAL tab_cntrl(length) ! tableau des parametres du run 620 character(len=*),parameter :: modname='dynredem1' 621 character*80 abort_message 622 c 623 INTEGER nb 624 SAVE nb 625 DATA nb / 0 / 626 character(len=12) :: start_file_type="earth" ! default start file type 627 628 if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then 629 write(lunout,*) trim(modname),' : Planeto-like start file' 630 start_file_type="planeto" 631 else 632 write(lunout,*) trim(modname),' : Earth-like start file' 633 endif 634 635 ierr = NF_OPEN(fichnom, NF_WRITE, nid) 636 IF (ierr .NE. NF_NOERR) THEN 637 PRINT*, "dynredem1: Pb. d ouverture "//trim(fichnom) 638 CALL abort_gcm("dynredem1", "", 1) 639 ENDIF 640 641 c Ecriture/extension de la coordonnee temps 642 643 nb = nb + 1 644 if (start_file_type.eq."earth") then 645 ierr = NF_INQ_VARID(nid, "temps", nvarid) 283 IMPLICIT NONE 284 include "dimensions.h" 285 include "paramet.h" 286 include "netcdf.inc" 287 include "comgeom.h" 288 include "iniprint.h" 289 !=============================================================================== 290 ! Arguments: 291 CHARACTER(LEN=*), INTENT(IN) :: fichnom !-- FILE NAME 292 REAL, INTENT(IN) :: time !-- TIME 293 REAL, INTENT(IN) :: vcov(iip1,jjm, llm) !-- V COVARIANT WIND 294 REAL, INTENT(IN) :: ucov(iip1,jjp1,llm) !-- U COVARIANT WIND 295 REAL, INTENT(IN) :: teta(iip1,jjp1,llm) !-- POTENTIAL TEMPERATURE 296 REAL, INTENT(INOUT) :: q(iip1,jjp1,llm,nqtot) !-- TRACERS 297 REAL, INTENT(IN) :: masse(iip1,jjp1,llm) !-- MASS PER CELL 298 REAL, INTENT(IN) :: ps(iip1,jjp1) !-- GROUND PRESSURE 299 !=============================================================================== 300 ! Local variables: 301 INTEGER :: l, iq, nid, vID, ierr, nid_trac, vID_trac 302 INTEGER,SAVE :: nb=0 303 INTEGER, PARAMETER :: length=100 304 REAL :: tab_cntrl(length) ! tableau des parametres du run 305 CHARACTER(LEN=256) :: var, dum 306 LOGICAL :: lread_inca 307 CHARACTER(LEN=80) :: abort_message 308 CHARACTER(len=12) :: start_file_type="earth" ! default start file type 309 310 ! fill dynredem_mod module variables 311 modname='dynredem1'; fil=fichnom 312 313 if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then 314 write(lunout,*) trim(modname),' : Planeto-like start file' 315 start_file_type="planeto" 316 else 317 write(lunout,*) trim(modname),' : Earth-like start file' 318 endif 319 320 CALL err(NF90_OPEN(fil,NF90_WRITE,nid),"open",fil) 321 322 !--- Write/extend time coordinate 323 nb = nb + 1 324 if (start_file_type.eq."earth") then 325 ierr = NF_INQ_VARID(nid, "temps", vID) 646 326 IF (ierr .NE. NF_NOERR) THEN 647 327 write(lunout,*) NF_STRERROR(ierr) … … 649 329 CALL abort_gcm(modname,abort_message,ierr) 650 330 ENDIF 651 652 ierr = NF_INQ_VARID(nid,"Time", nvarid)331 else 332 ierr = NF_INQ_VARID(nid,"Time", vID) 653 333 IF (ierr .NE. NF_NOERR) THEN 654 334 write(lunout,*) NF_STRERROR(ierr) … … 656 336 CALL abort_gcm(modname,abort_message,ierr) 657 337 ENDIF 658 endif ! of if (start_file_type.eq."earth") 659 call NF95_PUT_VAR(nid,nvarid,time,start=(/nb/)) 660 write(lunout,*) "dynredem1: Enregistrement pour ", nb, time 661 662 c 663 c Re-ecriture du tableau de controle, itaufin n'est plus defini quand 664 c on passe dans dynredem0 665 ierr = NF_INQ_VARID (nid, "controle", nvarid) 666 IF (ierr .NE. NF_NOERR) THEN 667 abort_message="dynredem1: Le champ <controle> est absent" 668 ierr = 1 669 CALL abort_gcm(modname,abort_message,ierr) 670 ENDIF 671 ierr = NF90_GET_VAR(nid, nvarid, tab_cntrl) 672 if (start_file_type=="earth") then 673 tab_cntrl(31) = REAL(itau_dyn + itaufin) 674 else 675 tab_cntrl(31) = 0 676 endif 677 call NF95_PUT_VAR(nid,nvarid,tab_cntrl) 678 679 c Ecriture des champs 680 c 681 ierr = NF_INQ_VARID(nid, "ucov", nvarid) 682 IF (ierr .NE. NF_NOERR) THEN 683 abort_message="Variable ucov n est pas definie" 684 ierr=1 685 CALL abort_gcm(modname,abort_message,ierr) 686 ENDIF 687 call NF95_PUT_VAR(nid,nvarid,ucov,start=(/1,1,1,nb/)) 688 689 ierr = NF_INQ_VARID(nid, "vcov", nvarid) 690 IF (ierr .NE. NF_NOERR) THEN 691 abort_message="Variable vcov n est pas definie" 692 ierr=1 693 CALL abort_gcm(modname,abort_message,ierr) 694 ENDIF 695 call NF95_PUT_VAR(nid,nvarid,vcov,start=(/1,1,1,nb/)) 696 697 ierr = NF_INQ_VARID(nid, "teta", nvarid) 698 IF (ierr .NE. NF_NOERR) THEN 699 abort_message="Variable teta n est pas definie" 700 ierr=1 701 CALL abort_gcm(modname,abort_message,ierr) 702 ENDIF 703 call NF95_PUT_VAR(nid,nvarid,teta,start=(/1,1,1,nb/)) 704 705 IF (type_trac == 'inca') THEN 706 ! Ajout Anne pour lecture valeurs traceurs dans un fichier start_trac.nc 707 ierr_file = NF_OPEN ("start_trac.nc", NF_NOWRITE,nid_trac) 708 IF (ierr_file .NE.NF_NOERR) THEN 709 write(lunout,*)'dynredem1: Pb d''ouverture du fichier', 710 & ' start_trac.nc' 711 write(lunout,*)' ierr = ', ierr_file 712 ENDIF 713 END IF 714 715 IF(nqtot.GE.1) THEN 716 do iq=1,nqtot 717 718 IF (type_trac /= 'inca') THEN 719 ierr = NF_INQ_VARID(nid, tname(iq), nvarid) 720 IF (ierr .NE. NF_NOERR) THEN 721 abort_message="Variable tname(iq) n est pas definie" 722 ierr=1 723 CALL abort_gcm(modname,abort_message,ierr) 724 ENDIF 725 call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq),start=(/1,1,1,nb/)) 726 ELSE ! type_trac = inca 727 ! lecture de la valeur du traceur dans start_trac.nc 728 IF (ierr_file .ne. 2) THEN 729 ierr = NF_INQ_VARID (nid_trac, tname(iq), nvarid_trac) 730 IF (ierr .NE. NF_NOERR) THEN 731 write(lunout,*) "dynredem1: ",trim(tname(iq)), 732 & " est absent de start_trac.nc" 733 ierr = NF_INQ_VARID(nid, tname(iq), nvarid) 734 IF (ierr .NE. NF_NOERR) THEN 735 abort_message="dynredem1: Variable "// 736 & trim(tname(iq))//" n est pas definie" 737 ierr=1 738 CALL abort_gcm(modname,abort_message,ierr) 739 ENDIF 740 call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq)) 741 742 ELSE 743 write(lunout,*) "dynredem1: ",trim(tname(iq)), 744 & " est present dans start_trac.nc" 745 ierr = NF90_GET_VAR(nid_trac, nvarid_trac, trac_tmp) 746 IF (ierr .NE. NF_NOERR) THEN 747 abort_message="dynredem1: Lecture echouee pour"// 748 & trim(tname(iq)) 749 ierr=1 750 CALL abort_gcm(modname,abort_message,ierr) 751 ENDIF 752 ierr = NF_INQ_VARID(nid, tname(iq), nvarid) 753 IF (ierr .NE. NF_NOERR) THEN 754 abort_message="dynredem1: Variable "// 755 & trim(tname(iq))//" n est pas definie" 756 ierr=1 757 CALL abort_gcm(modname,abort_message,ierr) 758 ENDIF 759 call NF95_PUT_VAR(nid, nvarid, trac_tmp) 760 761 ENDIF ! IF (ierr .NE. NF_NOERR) 762 ! fin lecture du traceur 763 ELSE ! si il n'y a pas de fichier start_trac.nc 764 ! print *, 'il n y a pas de fichier start_trac' 765 ierr = NF_INQ_VARID(nid, tname(iq), nvarid) 766 IF (ierr .NE. NF_NOERR) THEN 767 abort_message="dynredem1: Variable "// 768 & trim(tname(iq))//" n est pas definie" 769 ierr=1 770 CALL abort_gcm(modname,abort_message,ierr) 771 ENDIF 772 call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq), 773 & start=(/1,1,1,nb/)) 774 ENDIF ! (ierr_file .ne. 2) 775 END IF !type_trac 776 777 ENDDO 778 ENDIF 779 c 780 ierr = NF_INQ_VARID(nid, "masse", nvarid) 781 IF (ierr .NE. NF_NOERR) THEN 782 abort_message="dynredem1: Variable masse n est pas definie" 783 ierr=1 784 CALL abort_gcm(modname,abort_message,ierr) 785 ENDIF 786 call NF95_PUT_VAR(nid,nvarid,masse,start=(/1,1,1,nb/)) 787 c 788 ierr = NF_INQ_VARID(nid, "ps", nvarid) 789 IF (ierr .NE. NF_NOERR) THEN 790 abort_message="dynredem1: Variable ps n est pas definie" 791 ierr=1 792 CALL abort_gcm(modname,abort_message,ierr) 793 ENDIF 794 call NF95_PUT_VAR(nid,nvarid,ps,start=(/1,1,nb/)) 795 796 ierr = NF_CLOSE(nid) 797 c 798 RETURN 799 END 800 338 endif ! of if (start_file_type.eq."earth") 339 call NF95_PUT_VAR(nid,vID,time,start=(/nb/)) 340 WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time 341 342 !--- Rewrite control table (itaufin undefined in dynredem0) 343 var="controle" 344 CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var) 345 CALL err(NF90_GET_VAR(nid,vID,tab_cntrl),"get",var) 346 if (start_file_type=="earth") then 347 tab_cntrl(31) = REAL(itau_dyn + itaufin) 348 else 349 tab_cntrl(31) = 0 350 endif 351 CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var) 352 CALL err(NF90_PUT_VAR(nid,vID,tab_cntrl),"put",var) 353 354 !--- Save fields 355 CALL dynredem_write_u(nid,"ucov" ,ucov ,llm, nb) 356 CALL dynredem_write_v(nid,"vcov" ,vcov ,llm, nb) 357 CALL dynredem_write_u(nid,"teta" ,teta ,llm, nb) 358 CALL dynredem_write_u(nid,"masse",masse,llm, nb) 359 CALL dynredem_write_u(nid,"ps" ,ps ,1, nb) 360 361 !--- Tracers in file "start_trac.nc" (added by Anne) 362 lread_inca=.FALSE.; fil="start_trac.nc" 363 IF(type_trac=='inca') INQUIRE(FILE=fil,EXIST=lread_inca) 364 IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open") 365 366 !--- Save tracers 367 IF(nqtot.GE.1) THEN 368 DO iq=1,nqtot 369 var=tname(iq); ierr=-1 370 IF(lread_inca) THEN !--- Possibly read from "start_trac.nc" 371 fil="start_trac.nc" 372 ierr=NF90_INQ_VARID(nid_trac,var,vID_trac) 373 dum='inq'; IF(ierr==NF90_NoErr) dum='fnd' 374 WRITE(lunout,*)msg(dum,var) 375 376 377 IF(ierr==NF90_NoErr) CALL dynredem_read_u(nid_trac,var,q(:,:,:,iq),llm) 378 END IF ! of IF(lread_inca) 379 fil=fichnom 380 CALL dynredem_write_u(nid,var,q(:,:,:,iq),llm,nb) 381 END DO ! of DO iq=1,nqtot 382 ENDIF ! of IF(nqtot.GE.1) 383 384 CALL err(NF90_CLOSE(nid),"close") 385 fil="start_trac.nc" 386 IF(lread_inca) CALL err(NF90_CLOSE(nid_trac),"close") 387 388 END SUBROUTINE dynredem1 389 -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d_common/infotrac.F90 ¶
r1403 r1508 12 12 INTEGER, SAVE :: nbtr 13 13 14 ! CRisi: nb of father tracers (i.e. directly advected by air) 15 INTEGER, SAVE :: nqperes 16 14 17 ! Name variables 15 18 CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics … … 22 25 ! dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code. 23 26 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique 27 28 ! CRisi: arrays for sons 29 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqfils 30 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: nqdesc ! number of sons + all gran-sons over all generations 31 INTEGER, SAVE :: nqdesc_tot 32 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqfils 33 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iqpere 24 34 25 35 ! conv_flg(it)=0 : convection desactivated for tracer number it … … 30 40 CHARACTER(len=4),SAVE :: type_trac 31 41 CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym 32 42 43 ! CRisi: specific stuff for isotopes 44 LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso 45 INTEGER :: niso_possibles 46 PARAMETER ( niso_possibles=5) ! 5 possible water isotopes 47 real, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal 48 LOGICAL, DIMENSION(niso_possibles),SAVE :: use_iso 49 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase) 50 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot 51 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot 52 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numéro de la zone de tracage en fn de nqtot 53 INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numéro de la zone de tracage en fn de nqtot 54 INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles 55 INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: index_trac ! numéro ixt en fn izone, indnum entre 1 et niso 56 INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso 57 33 58 CONTAINS 34 59 35 60 SUBROUTINE infotrac_init 36 USE control_mod 61 USE control_mod, ONLY: planet_type, config_inca 37 62 #ifdef REPROBUS 38 63 USE CHEM_REP, ONLY : Init_chem_rep_trac … … 63 88 64 89 CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_0 ! tracer short name 90 CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi 65 91 CHARACTER(len=3), DIMENSION(30) :: descrq 66 92 CHARACTER(len=1), DIMENSION(3) :: txts … … 70 96 INTEGER :: nqtrue ! number of tracers read from tracer.def, without higer order of moment 71 97 INTEGER :: iq, new_iq, iiq, jq, ierr, ierr2, ierr3 98 INTEGER :: ifils,ipere,generation ! CRisi 99 LOGICAL :: continu,nouveau_traceurdef 100 INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def 101 CHARACTER(len=15) :: tchaine 72 102 73 103 character(len=80) :: line ! to store a line of text … … 138 168 WRITE(lunout,*) trim(modname),': Open traceur.def : ok' 139 169 READ(90,*) nqtrue 170 WRITE(lunout,*) trim(modname),' nqtrue=',nqtrue 140 171 ELSE 141 172 WRITE(lunout,*) trim(modname),': Problem in opening traceur.def' … … 146 177 nbtr=nqtrue-2 147 178 ELSE ! type_trac=inca 148 ! nbtr has been read from INCA by init_cont_lmdz() in gcm.F 149 nqtrue=nbtr+2 150 END IF 179 ! The traceur.def file is used to define the number "nqo" of water phases 180 ! present in the simulation. Default : nqo = 2. 181 OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr) 182 IF(ierr.EQ.0) THEN 183 WRITE(lunout,*) trim(modname),': Open traceur.def : ok' 184 READ(90,*) nqo 185 ELSE 186 WRITE(lunout,*) trim(modname),': Using default value for nqo' 187 nqo=2 188 ENDIF 189 IF (nqo /= 2 .AND. nqo /= 3 ) THEN 190 WRITE(lunout,*) trim(modname),': nqo=',nqo, ' is not allowded. Only 2 or 3 water phases allowed' 191 CALL abort_gcm('infotrac_init','Bad number of water phases',1) 192 END IF 193 ! nbtr has been read from INCA by init_const_lmdz() in gcm.F 194 #ifdef INCA 195 CALL Init_chem_inca_trac(nbtr) 196 #endif 197 nqtrue=nbtr+nqo 198 END IF ! type_trac 151 199 152 200 IF (nqtrue < 2) THEN … … 155 203 END IF 156 204 205 !jyg< 157 206 ! Transfert number of tracers to Reprobus 158 IF (type_trac == 'repr') THEN 159 #ifdef REPROBUS 160 CALL Init_chem_rep_trac(nbtr) 161 #endif 162 END IF 207 !! IF (type_trac == 'repr') THEN 208 !!#ifdef REPROBUS 209 !! CALL Init_chem_rep_trac(nbtr) 210 !!#endif 211 !! END IF 212 !>jyg 163 213 164 214 ELSE ! not Earth … … 175 225 ! in the dynamics than in the physics 176 226 nbtr=nqtrue 227 nqo=0 177 228 178 229 ENDIF ! planet_type … … 180 231 ! Allocate variables depending on nqtrue and nbtr 181 232 ! 182 ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue)) 183 ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr)) 184 conv_flg(:) = 1 ! convection activated for all tracers 185 pbl_flg(:) = 1 ! boundary layer activated for all tracers 233 ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue)) 234 ! 235 !jyg< 236 !! ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr)) 237 !! conv_flg(:) = 1 ! convection activated for all tracers 238 !! pbl_flg(:) = 1 ! boundary layer activated for all tracers 239 !>jyg 186 240 187 241 !----------------------------------------------------------------------- … … 216 270 ! Continue to read tracer.def 217 271 DO iq=1,nqtrue 218 READ(90,*) hadv(iq),vadv(iq),tnom_0(iq) 219 END DO 272 273 write(*,*) 'infotrac 237: iq=',iq 274 ! CRisi: ajout du nom du fluide transporteur 275 ! mais rester retro compatible 276 READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine 277 write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq) 278 write(lunout,*) 'tchaine=',trim(tchaine) 279 ! write(*,*) 'infotrac 238: IOstatus=',IOstatus 280 if (IOstatus.ne.0) then 281 CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1) 282 endif 283 ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un 284 ! espace ou pas au milieu de la chaine. 285 continu=.true. 286 nouveau_traceurdef=.false. 287 iiq=1 288 do while (continu) 289 if (tchaine(iiq:iiq).eq.' ') then 290 nouveau_traceurdef=.true. 291 continu=.false. 292 else if (iiq.lt.LEN_TRIM(tchaine)) then 293 iiq=iiq+1 294 else 295 continu=.false. 296 endif 297 enddo 298 write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef 299 if (nouveau_traceurdef) then 300 write(lunout,*) 'C''est la nouvelle version de traceur.def' 301 tnom_0(iq)=tchaine(1:iiq-1) 302 tnom_transp(iq)=tchaine(iiq+1:15) 303 else 304 write(lunout,*) 'C''est l''ancienne version de traceur.def' 305 write(lunout,*) 'On suppose que les traceurs sont tous d''air' 306 tnom_0(iq)=tchaine 307 tnom_transp(iq) = 'air' 308 endif 309 write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>' 310 write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>' 311 312 END DO ! DO iq=1,nqtrue 220 313 CLOSE(90) 314 221 315 ELSE ! Without tracer.def, set default values (for Earth!) 222 316 if ((nqtrue==4).and.(planet_type=="earth")) then … … 224 318 vadv(1) = 14 225 319 tnom_0(1) = 'H2Ov' 320 tnom_transp(1) = 'air' 226 321 hadv(2) = 10 227 322 vadv(2) = 10 228 323 tnom_0(2) = 'H2Ol' 324 tnom_transp(2) = 'air' 229 325 hadv(3) = 10 230 326 vadv(3) = 10 231 327 tnom_0(3) = 'RN' 328 tnom_transp(3) = 'air' 232 329 hadv(4) = 10 233 330 vadv(4) = 10 234 331 tnom_0(4) = 'PB' 332 tnom_transp(4) = 'air' 235 333 else 236 334 ! Error message, we need a traceur.def file … … 240 338 CALL abort_gcm('infotrac_init','Need a traceur.def file!',1) 241 339 endif ! of if (nqtrue==4) 242 END IF 340 END IF ! of IF(ierr.EQ.0) 243 341 244 !CR: nombre de traceurs de l eau245 if (tnom_0(3) == 'H2Oi') then246 nqo=3247 else248 nqo=2249 endif250 251 342 WRITE(lunout,*) trim(modname),': Valeur de traceur.def :' 252 343 WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue 253 344 DO iq=1,nqtrue 254 WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq) 345 WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq) 255 346 END DO 256 347 257 ELSE ! type_trac=inca : config_inca='aero' ou 'chem' 348 !CR: nombre de traceurs de l eau 349 if (tnom_0(3) == 'H2Oi') then 350 nqo=3 351 else 352 nqo=2 353 endif 354 ! For Earth, water vapour & liquid tracers are not in the physics 355 nbtr=nqtrue-nqo 356 ENDIF ! (type_trac == 'lmdz' .OR. type_trac == 'repr') 357 !jyg< 358 ! 359 ! Transfert number of tracers to Reprobus 360 IF (type_trac == 'repr') THEN 361 #ifdef REPROBUS 362 CALL Init_chem_rep_trac(nbtr) 363 #endif 364 END IF 365 ! 366 ! Allocate variables depending on nbtr 367 ! 368 ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr)) 369 conv_flg(:) = 1 ! convection activated for all tracers 370 pbl_flg(:) = 1 ! boundary layer activated for all tracers 371 ! 372 !! ELSE ! type_trac=inca : config_inca='aero' ou 'chem' 373 ! 374 IF (type_trac == 'inca') THEN ! config_inca='aero' ou 'chem' 375 !>jyg 258 376 ! le module de chimie fournit les noms des traceurs 259 377 ! et les schemas d'advection associes. … … 268 386 #endif 269 387 tnom_0(1)='H2Ov' 388 tnom_transp(1) = 'air' 270 389 tnom_0(2)='H2Ol' 271 272 DO iq =3,nqtrue 273 tnom_0(iq)=solsym(iq-2) 390 tnom_transp(2) = 'air' 391 IF (nqo == 3) then 392 tnom_0(3)='H2Oi' !! jyg 393 tnom_transp(3) = 'air' 394 endif 395 396 !jyg< 397 DO iq = nqo+1, nqtrue 398 tnom_0(iq)=solsym(iq-nqo) 399 tnom_transp(iq) = 'air' 274 400 END DO 275 nqo = 2 276 277 END IF ! type_trac 401 !! DO iq =3,nqtrue 402 !! tnom_0(iq)=solsym(iq-2) 403 !! END DO 404 !! nqo = 2 405 !>jyg 406 407 END IF ! (type_trac == 'inca') 278 408 279 409 ELSE ! not Earth 410 ! Other planets (for now); we have the same number of tracers 411 ! in the dynamics than in the physics 412 nbtr=nqtrue 413 ! NB: Reading a traceur.def with isotopes remains to be done... 280 414 IF(ierr.EQ.0) THEN 281 415 ! Continue to read tracer.def … … 296 430 endif 297 431 endif ! of if(ierr2.ne.0) 432 tnom_transp(iq)='air' ! no isotopes... for now... 298 433 END DO ! of DO iq=1,nqtrue 299 434 CLOSE(90) … … 302 437 vadv(1) = 10 303 438 tnom_0(1) = 'dummy' 439 tnom_transp(1)='air' 304 440 END IF 305 441 … … 307 443 WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue 308 444 DO iq=1,nqtrue 309 WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq) 445 WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq) 310 446 END DO 311 447 … … 437 573 438 574 !----------------------------------------------------------------------- 575 ! 576 ! 5) Determine father/son relations for isotopes and carrying fluid 577 ! 578 !----------------------------------------------------------------------- 579 580 ! CRisi: quels sont les traceurs fils et les traceurs pères. 581 ! initialiser tous les tableaux d'indices liés aux traceurs familiaux 582 ! + vérifier que tous les pères sont écrits en premières positions 583 ALLOCATE(nqfils(nqtot),nqdesc(nqtot)) 584 ALLOCATE(iqfils(nqtot,nqtot)) 585 ALLOCATE(iqpere(nqtot)) 586 nqperes=0 587 nqfils(:)=0 588 nqdesc(:)=0 589 iqfils(:,:)=0 590 iqpere(:)=0 591 nqdesc_tot=0 592 DO iq=1,nqtot 593 if (tnom_transp(iq) == 'air') then 594 ! ceci est un traceur père 595 WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere' 596 nqperes=nqperes+1 597 iqpere(iq)=0 598 else !if (tnom_transp(iq) == 'air') then 599 ! ceci est un fils. Qui est son père? 600 WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un fils' 601 continu=.true. 602 ipere=1 603 do while (continu) 604 if (tnom_transp(iq) == tnom_0(ipere)) then 605 ! Son père est ipere 606 WRITE(lunout,*) 'Le traceur',iq,'appele ', & 607 & trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere)) 608 nqfils(ipere)=nqfils(ipere)+1 609 iqfils(nqfils(ipere),ipere)=iq 610 iqpere(iq)=ipere 611 continu=.false. 612 else !if (tnom_transp(iq) == tnom_0(ipere)) then 613 ipere=ipere+1 614 if (ipere.gt.nqtot) then 615 WRITE(lunout,*) 'Le traceur',iq,'appele ', & 616 & trim(tnom_0(iq)),', est orpelin.' 617 CALL abort_gcm('infotrac_init','Un traceur est orphelin',1) 618 endif !if (ipere.gt.nqtot) then 619 endif !if (tnom_transp(iq) == tnom_0(ipere)) then 620 enddo !do while (continu) 621 endif !if (tnom_transp(iq) == 'air') then 622 enddo !DO iq=1,nqtot 623 WRITE(lunout,*) 'infotrac: nqperes=',nqperes 624 WRITE(lunout,*) 'nqfils=',nqfils 625 WRITE(lunout,*) 'iqpere=',iqpere 626 WRITE(lunout,*) 'iqfils=',iqfils 627 628 ! Calculer le nombre de descendants à partir de iqfils et de nbfils 629 DO iq=1,nqtot 630 generation=0 631 continu=.true. 632 ifils=iq 633 do while (continu) 634 ipere=iqpere(ifils) 635 if (ipere.gt.0) then 636 nqdesc(ipere)=nqdesc(ipere)+1 637 nqdesc_tot=nqdesc_tot+1 638 iqfils(nqdesc(ipere),ipere)=iq 639 ifils=ipere 640 generation=generation+1 641 else !if (ipere.gt.0) then 642 continu=.false. 643 endif !if (ipere.gt.0) then 644 enddo !do while (continu) 645 WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation 646 enddo !DO iq=1,nqtot 647 WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc 648 WRITE(lunout,*) 'iqfils=',iqfils 649 WRITE(lunout,*) 'nqdesc_tot=',nqdesc_tot 650 651 ! Interdire autres schémas que 10 pour les traceurs fils, et autres schémas 652 ! que 10 et 14 si des pères ont des fils 653 do iq=1,nqtot 654 if (iqpere(iq).gt.0) then 655 ! ce traceur a un père qui n'est pas l'air 656 ! Seul le schéma 10 est autorisé 657 if (iadv(iq)/=10) then 658 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons' 659 CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1) 660 endif 661 ! Le traceur père ne peut être advecté que par schéma 10 ou 14: 662 IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN 663 WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for fathers' 664 CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1) 665 endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN 666 endif !if (iqpere(iq).gt.0) the 667 enddo !do iq=1,nqtot 668 669 670 ! detecter quels sont les traceurs isotopiques parmi des traceurs 671 call infotrac_isoinit(tnom_0,nqtrue) 672 673 !----------------------------------------------------------------------- 439 674 ! Finalize : 440 675 ! 441 DEALLOCATE(tnom_0, hadv, vadv )676 DEALLOCATE(tnom_0, hadv, vadv,tnom_transp) 442 677 443 678 444 679 END SUBROUTINE infotrac_init 680 681 !----------------------------------------------------------------------- 682 683 SUBROUTINE infotrac_isoinit(tnom_0,nqtrue) 684 685 #ifdef CPP_IOIPSL 686 use IOIPSL 687 #else 688 ! if not using IOIPSL, we still need to use (a local version of) getin 689 use ioipsl_getincom 690 #endif 691 implicit none 692 693 ! inputs 694 INTEGER nqtrue 695 CHARACTER(len=15) tnom_0(nqtrue) 696 697 ! locals 698 CHARACTER(len=3), DIMENSION(niso_possibles) :: tnom_iso 699 INTEGER, ALLOCATABLE,DIMENSION(:,:) :: nb_iso,nb_traciso 700 INTEGER, ALLOCATABLE,DIMENSION(:) :: nb_isoind 701 INTEGER :: ntraceurs_zone_prec,iq,phase,ixt,iiso,izone 702 CHARACTER(len=19) :: tnom_trac 703 INCLUDE "iniprint.h" 704 705 tnom_iso=(/'eau','HDO','O18','O17','HTO'/) 706 707 ALLOCATE(nb_iso(niso_possibles,nqo)) 708 ALLOCATE(nb_isoind(nqo)) 709 ALLOCATE(nb_traciso(niso_possibles,nqo)) 710 ALLOCATE(iso_num(nqtot)) 711 ALLOCATE(iso_indnum(nqtot)) 712 ALLOCATE(zone_num(nqtot)) 713 ALLOCATE(phase_num(nqtot)) 714 715 iso_num(:)=0 716 iso_indnum(:)=0 717 zone_num(:)=0 718 phase_num(:)=0 719 indnum_fn_num(:)=0 720 use_iso(:)=.false. 721 nb_iso(:,:)=0 722 nb_isoind(:)=0 723 nb_traciso(:,:)=0 724 niso=0 725 ntraceurs_zone=0 726 ntraceurs_zone_prec=0 727 ntraciso=0 728 729 do iq=nqo+1,nqtot 730 write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq) 731 do phase=1,nqo 732 do ixt= 1,niso_possibles 733 tnom_trac=trim(tnom_0(phase))//'_' 734 tnom_trac=trim(tnom_trac)//trim(tnom_iso(ixt)) 735 write(*,*) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac 736 IF (tnom_0(iq) == tnom_trac) then 737 write(lunout,*) 'Ce traceur est un isotope' 738 nb_iso(ixt,phase)=nb_iso(ixt,phase)+1 739 nb_isoind(phase)=nb_isoind(phase)+1 740 iso_num(iq)=ixt 741 iso_indnum(iq)=nb_isoind(phase) 742 indnum_fn_num(ixt)=iso_indnum(iq) 743 phase_num(iq)=phase 744 write(lunout,*) 'iso_num(iq)=',iso_num(iq) 745 write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq) 746 write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt) 747 write(lunout,*) 'phase_num(iq)=',phase_num(iq) 748 goto 20 749 else if (iqpere(iq).gt.0) then 750 if (tnom_0(iqpere(iq)) == tnom_trac) then 751 write(lunout,*) 'Ce traceur est le fils d''un isotope' 752 ! c'est un traceur d'isotope 753 nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1 754 iso_num(iq)=ixt 755 iso_indnum(iq)=indnum_fn_num(ixt) 756 zone_num(iq)=nb_traciso(ixt,phase) 757 phase_num(iq)=phase 758 write(lunout,*) 'iso_num(iq)=',iso_num(iq) 759 write(lunout,*) 'phase_num(iq)=',phase_num(iq) 760 write(lunout,*) 'zone_num(iq)=',zone_num(iq) 761 goto 20 762 endif !if (tnom_0(iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then 763 endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then 764 enddo !do ixt= niso_possibles 765 enddo !do phase=1,nqo 766 20 continue 767 enddo !do iq=1,nqtot 768 769 write(lunout,*) 'iso_num=',iso_num 770 write(lunout,*) 'iso_indnum=',iso_indnum 771 write(lunout,*) 'zone_num=',zone_num 772 write(lunout,*) 'phase_num=',phase_num 773 write(lunout,*) 'indnum_fn_num=',indnum_fn_num 774 775 do ixt= 1,niso_possibles 776 777 if (nqo.gt.0) then ! Ehouarn: because tests below only valid if nqo>=1, 778 779 if (nb_iso(ixt,1).eq.1) then 780 ! on vérifie que toutes les phases ont le même nombre de 781 ! traceurs 782 do phase=2,nqo 783 if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then 784 ! write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase) 785 CALL abort_gcm('infotrac_init','Phases must have same number of isotopes',1) 786 endif 787 enddo !do phase=2,nqo 788 789 niso=niso+1 790 use_iso(ixt)=.true. 791 ntraceurs_zone=nb_traciso(ixt,1) 792 793 ! on vérifie que toutes les phases ont le même nombre de 794 ! traceurs 795 do phase=2,nqo 796 if (nb_traciso(ixt,phase).ne.ntraceurs_zone) then 797 write(lunout,*) 'ixt,phase,nb_traciso=',ixt,phase,nb_traciso(ixt,phase) 798 write(lunout,*) 'ntraceurs_zone=',ntraceurs_zone 799 CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1) 800 endif 801 enddo !do phase=2,nqo 802 ! on vérifie que tous les isotopes ont le même nombre de 803 ! traceurs 804 if (ntraceurs_zone_prec.gt.0) then 805 if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 806 ntraceurs_zone_prec=ntraceurs_zone 807 else !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 808 write(*,*) 'ntraceurs_zone_prec,ntraceurs_zone=',ntraceurs_zone_prec,ntraceurs_zone 809 CALL abort_gcm('infotrac_init', & 810 &'Isotope tracers are not well defined in traceur.def',1) 811 endif !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then 812 endif !if (ntraceurs_zone_prec.gt.0) then 813 814 else if (nb_iso(ixt,1).ne.0) then 815 WRITE(lunout,*) 'nqo,ixt=',nqo,ixt 816 WRITE(lunout,*) 'nb_iso(ixt,1)=',nb_iso(ixt,1) 817 CALL abort_gcm('infotrac_init','Isotopes are not well defined in traceur.def',1) 818 endif !if (nb_iso(ixt,1).eq.1) then 819 820 endif ! of if (nqo.gt.0) 821 822 enddo ! do ixt= niso_possibles 823 824 ! dimensions isotopique: 825 ntraciso=niso*(ntraceurs_zone+1) 826 WRITE(lunout,*) 'niso=',niso 827 WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso 828 829 ! flags isotopiques: 830 if (niso.gt.0) then 831 ok_isotopes=.true. 832 else 833 ok_isotopes=.false. 834 endif 835 WRITE(lunout,*) 'ok_isotopes=',ok_isotopes 836 837 if (ok_isotopes) then 838 ok_iso_verif=.false. 839 call getin('ok_iso_verif',ok_iso_verif) 840 ok_init_iso=.false. 841 call getin('ok_init_iso',ok_init_iso) 842 tnat=(/1.0,155.76e-6,2005.2e-6,0.004/100.,0.0/) 843 alpha_ideal=(/1.0,1.01,1.006,1.003,1.0/) 844 endif !if (ok_isotopes) then 845 WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif 846 WRITE(lunout,*) 'ok_init_iso=',ok_init_iso 847 848 if (ntraceurs_zone.gt.0) then 849 ok_isotrac=.true. 850 else 851 ok_isotrac=.false. 852 endif 853 WRITE(lunout,*) 'ok_isotrac=',ok_isotrac 854 855 ! remplissage du tableau iqiso(ntraciso,phase) 856 ALLOCATE(iqiso(ntraciso,nqo)) 857 iqiso(:,:)=0 858 do iq=1,nqtot 859 if (iso_num(iq).gt.0) then 860 ixt=iso_indnum(iq)+zone_num(iq)*niso 861 iqiso(ixt,phase_num(iq))=iq 862 endif 863 enddo 864 ! WRITE(lunout,*) 'iqiso=',iqiso 865 866 ! replissage du tableau index_trac(ntraceurs_zone,niso) 867 ALLOCATE(index_trac(ntraceurs_zone,niso)) 868 if (ok_isotrac) then 869 do iiso=1,niso 870 do izone=1,ntraceurs_zone 871 index_trac(izone,iiso)=iiso+izone*niso 872 enddo 873 enddo 874 else !if (ok_isotrac) then 875 index_trac(:,:)=0.0 876 endif !if (ok_isotrac) then 877 write(lunout,*) 'index_trac=',index_trac 878 879 ! Finalize : 880 DEALLOCATE(nb_iso) 881 882 END SUBROUTINE infotrac_isoinit 883 884 !----------------------------------------------------------------------- 445 885 446 886 ! Ehouarn: routine iniadvtrac => from Mars/generic; does essentially the -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d_common/iniacademic.F90 ¶
r1422 r1508 5 5 6 6 USE filtreg_mod, ONLY: inifilr 7 USE infotrac, ONLY : nqtot 7 USE infotrac, ONLY : nqtot, ok_isotopes, iso_num, zone_num, & 8 iqpere, tnat, alpha_ideal, iso_indnum, & 9 phase_num, iqiso, ok_iso_verif 8 10 USE control_mod, ONLY: day_step,planet_type 9 11 #ifdef CPP_IOIPSL … … 262 264 if (i == 2) q(:,:,i)=1.e-15 263 265 if (i.gt.2) q(:,:,i)=0. 266 267 ! CRisi: init des isotopes 268 ! distill de Rayleigh très simplifiée 269 if (ok_isotopes) then 270 if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then 271 q(:,:,i)=q(:,:,iqpere(i)) & 272 & *tnat(iso_num(i)) & 273 & *(q(:,:,iqpere(i))/30.e-3) & 274 & **(alpha_ideal(iso_num(i))-1) 275 endif 276 if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then 277 q(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i))) 278 endif 279 endif !if (ok_isotopes) then 280 264 281 enddo 265 282 else 266 283 q(:,:,:)=0 267 284 endif ! of if (planet_type=="earth") 285 286 if (ok_iso_verif) then 287 ! Ehouarn: this will onyly work in serial mode 288 ! call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc') 289 endif !if (ok_iso_verif) then 268 290 269 291 ! add random perturbation to temperature -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3d_common/pression.F90 ¶
r1506 r1508 1 SUBROUTINE pression( ngrid, ap, bp, ps, p ) 1 2 ! 2 ! $Header$ 3 ! 4 SUBROUTINE pression( ngrid, ap, bp, ps, p ) 5 c 3 !------------------------------------------------------------------------------- 4 ! Authors: P. Le Van , Fr.Hourdin 5 !------------------------------------------------------------------------------- 6 ! Purpose: Compute pressure p(l) at different levels from l = 1 (ground level) 7 ! to l = llm +1. Those levels correspond to the llm layers interfaces, 8 ! with p(ij,llm+1) = 0. and p(ij,1) = ps(ij) . 9 !------------------------------------------------------------------------------- 10 IMPLICIT NONE 11 include "dimensions.h" 12 include "paramet.h" 13 !=============================================================================== 14 ! Arguments: 15 INTEGER, INTENT(IN) :: ngrid !--- NUMBER OF GRID POINTS 16 REAL, INTENT(IN) :: ap(llmp1), bp(llmp1) !--- HYBRID COEFFICIENTS 17 REAL, INTENT(IN) :: ps(ngrid) !--- SURFACE PRESSURE 18 REAL, INTENT(OUT) :: p(ngrid,llmp1) !--- 3D PRESSURE FIELD 19 !=============================================================================== 20 ! Local variables: 21 INTEGER :: l 22 !=============================================================================== 23 DO l=1,llmp1; p(:,l) = ap(l) + bp(l) * ps(:); END DO 6 24 7 c Auteurs : P. Le Van , Fr.Hourdin . 25 END SUBROUTINE pression 8 26 9 c ************************************************************************ 10 c Calcule la pression p(l) aux differents niveaux l = 1 ( niveau du 11 c sol) a l = llm +1 ,ces niveaux correspondant aux interfaces des (llm) 12 c couches , avec p(ij,llm +1) = 0. et p(ij,1) = ps(ij) . 13 c ************************************************************************ 14 c 15 IMPLICIT NONE 16 c 17 #include "dimensions.h" 18 #include "paramet.h" 19 c 20 INTEGER ngrid 21 INTEGER l,ij 22 23 REAL ap( llmp1 ), bp( llmp1 ), ps( ngrid ), p( ngrid,llmp1 ) 24 25 DO l = 1, llmp1 26 DO ij = 1, ngrid 27 p(ij,l) = ap(l) + bp(l) * ps(ij) 28 ENDDO 29 ENDDO 30 31 RETURN 32 END 27 -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F90 ¶
r1507 r1508 2 2 ! $Id: dynredem_p.F 1635 2012-07-12 11:37:16Z lguez $ 3 3 ! 4 c 5 SUBROUTINE dynredem0_p(fichnom,iday_end,phis) 4 SUBROUTINE dynredem0_p(fichnom,iday_end,phis) 6 5 #ifdef CPP_IOIPSL 7 6 USE IOIPSL 8 7 #endif 9 USE parallel_lmdz 10 USE infotrac 11 use netcdf95, only: NF95_PUT_VAR 12 use control_mod, only : planet_type 13 USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,nivsig,nivsigs, 14 . presnivs,pseudoalt 15 USE comconst_mod, ONLY: daysec,dtvr,rad,omeg,g,cpp,kappa,pi 16 USE logic_mod, ONLY: fxyhypb,ysinus 17 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, 18 . taux,tauy 19 USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn,itaufin, 20 . start_time,hour_ini 21 USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 22 23 IMPLICIT NONE 24 c======================================================================= 25 c Ecriture du fichier de redemarrage sous format NetCDF (initialisation) 26 c======================================================================= 27 c Declarations: 28 c ------------- 29 #include "dimensions.h" 30 #include "paramet.h" 31 #include "comgeom2.h" 32 #include "netcdf.inc" 33 #include "iniprint.h" 34 35 c Arguments: 36 c ---------- 37 INTEGER iday_end 38 REAL phis(iip1, jjp1) 39 CHARACTER*(*) fichnom 40 41 c Local: 42 c ------ 43 INTEGER iq,l 44 INTEGER length 45 PARAMETER (length = 100) 46 REAL tab_cntrl(length) ! tableau des parametres du run 47 INTEGER ierr 48 character*20 modname 49 character*80 abort_message 50 51 c Variables locales pour NetCDF: 52 c 53 INTEGER dims2(2), dims3(3), dims4(4) 54 INTEGER idim_index 55 INTEGER idim_rlonu, idim_rlonv, idim_rlatu, idim_rlatv 56 INTEGER idim_s, idim_sig 57 INTEGER idim_tim 58 INTEGER nid,nvarid 59 60 REAL zan0,zjulian,hours 61 INTEGER yyears0,jjour0, mmois0 62 character*30 unites 63 64 character(len=12) :: start_file_type="earth" ! default start file type 65 INTEGER idecal 66 67 c----------------------------------------------------------------------- 68 if (mpi_rank==0) then 69 70 modname='dynredem0_p' 8 USE parallel_lmdz, ONLY: mpi_rank 9 USE infotrac, ONLY: nqtot, tname, ttext 10 USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL, & 11 NF90_CLOSE, NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER 12 USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil 13 use netcdf95, only: NF95_PUT_VAR 14 use control_mod, only : planet_type 15 USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt,pa,preff, & 16 nivsig,nivsigs 17 USE comconst_mod, ONLY: daysec,dtvr,rad,omeg,g,cpp,kappa,pi 18 USE logic_mod, ONLY: fxyhypb,ysinus 19 USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, & 20 taux,tauy 21 USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn,itaufin, & 22 start_time,hour_ini 23 USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0 24 25 IMPLICIT NONE 26 !======================================================================= 27 ! Writting the NetCDF restart file (initialisation) 28 !======================================================================= 29 ! Declarations: 30 ! ------------- 31 include "dimensions.h" 32 include "paramet.h" 33 include "comgeom2.h" 34 include "netcdf.inc" 35 include "iniprint.h" 36 37 !=============================================================================== 38 ! Arguments: 39 CHARACTER(LEN=*), INTENT(IN) :: fichnom !--- FILE NAME 40 INTEGER, INTENT(IN) :: iday_end !--- 41 REAL, INTENT(IN) :: phis(iip1, jjp1) !--- GROUND GEOPOTENTIAL 42 !=============================================================================== 43 ! Local variables: 44 INTEGER :: iq,l 45 INTEGER, PARAMETER :: length=100 46 REAL :: tab_cntrl(length) ! run parameters 47 INTEGER :: ierr 48 CHARACTER(LEN=80) :: abort_message 49 50 ! For NetCDF: 51 CHARACTER(LEN=30) :: unites 52 INTEGER :: indexID 53 INTEGER :: rlonuID, rlonvID, rlatuID, rlatvID 54 INTEGER :: sID, sigID, nID, vID, timID 55 INTEGER :: yyears0, jjour0, mmois0 56 REAL :: zan0, zjulian, hours 57 58 CHARACTER(len=12) :: start_file_type="earth" ! default start file type 59 INTEGER :: idecal 60 61 !=============================================================================== 62 if (mpi_rank==0) then ! only the master reads input file 63 ! fill dynredem_mod module variables 64 modname='dynredem0_p'; fil=fichnom 71 65 72 66 #ifdef CPP_IOIPSL 73 74 67 call ymds2ju(annee_ref, 1, iday_end, 0.0, zjulian) 68 call ju2ymds(zjulian, yyears0, mmois0, jjour0, hours) 75 69 #else 76 70 ! set yyears0, mmois0, jjour0 to 0,1,1 (hours is not used) 77 78 79 71 yyears0=0 72 mmois0=1 73 jjour0=1 80 74 #endif 81 75 82 !!! AS: idecal is a hack to be able to read planeto starts... 83 !!! .... while keeping everything OK for LMDZ EARTH 84 if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then 85 write(lunout,*) trim(modname),' : Planeto-like start file' 86 start_file_type="planeto" 87 idecal = 4 88 else 89 write(lunout,*) trim(modname),' : Earth-like start file' 90 idecal = 5 91 endif 92 93 DO l=1,length 94 tab_cntrl(l) = 0. 95 ENDDO 96 tab_cntrl(1) = REAL(iim) 97 tab_cntrl(2) = REAL(jjm) 98 tab_cntrl(3) = REAL(llm) 99 if (start_file_type.eq."earth") then 100 tab_cntrl(4)=REAL(day_ref) 101 else 102 !tab_cntrl(4)=REAL(day_end) 103 tab_cntrl(4)=REAL(iday_end) 104 endif 105 tab_cntrl(5) = REAL(annee_ref) 106 tab_cntrl(idecal+1) = rad 107 tab_cntrl(idecal+2) = omeg 108 tab_cntrl(idecal+3) = g 109 tab_cntrl(idecal+4) = cpp 110 tab_cntrl(idecal+5) = kappa 111 tab_cntrl(idecal+6) = daysec 112 tab_cntrl(idecal+7) = dtvr 113 tab_cntrl(idecal+8) = etot0 114 tab_cntrl(idecal+9) = ptot0 115 tab_cntrl(idecal+10) = ztot0 116 tab_cntrl(idecal+11) = stot0 117 tab_cntrl(idecal+12) = ang0 118 tab_cntrl(idecal+13) = pa 119 tab_cntrl(idecal+14) = preff 120 c 121 c ..... parametres pour le zoom ...... 122 123 tab_cntrl(idecal+15) = clon 124 tab_cntrl(idecal+16) = clat 125 tab_cntrl(idecal+17) = grossismx 126 tab_cntrl(idecal+18) = grossismy 127 c 128 IF ( fxyhypb ) THEN 129 tab_cntrl(idecal+19) = 1. 130 tab_cntrl(idecal+20) = dzoomx 131 tab_cntrl(idecal+21) = dzoomy 132 tab_cntrl(idecal+22) = 0. 133 tab_cntrl(idecal+23) = taux 134 tab_cntrl(idecal+24) = tauy 135 ELSE 136 tab_cntrl(idecal+19) = 0. 137 tab_cntrl(idecal+20) = dzoomx 138 tab_cntrl(idecal+21) = dzoomy 139 tab_cntrl(idecal+22) = 0. 140 tab_cntrl(idecal+23) = 0. 141 tab_cntrl(idecal+24) = 0. 142 IF( ysinus ) tab_cntrl(idecal+22) = 1. 143 ENDIF 144 145 if (start_file_type.eq."earth") then 146 tab_cntrl(idecal+25) = REAL(iday_end) 147 tab_cntrl(idecal+26) = REAL(itau_dyn + itaufin) 148 c start_time: start_time of simulation (not necessarily 0.) 149 tab_cntrl(idecal+27) = start_time 150 endif 151 152 if (planet_type=="mars") then ! For Mars only 153 tab_cntrl(29)=hour_ini 154 endif 155 c 156 c ......................................................... 157 c 158 c Creation du fichier: 159 c 160 ierr = NF_CREATE(fichnom, NF_CLOBBER, nid) 161 IF (ierr.NE.NF_NOERR) THEN 162 WRITE(6,*)" Pb d ouverture du fichier "//fichnom 163 WRITE(6,*)' ierr = ', ierr 164 CALL ABORT 165 ENDIF 166 c 167 c Preciser quelques attributs globaux: 168 c 169 ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 27, 170 . "Fichier demmarage dynamique") 171 c 172 c Definir les dimensions du fichiers: 173 c 174 if (start_file_type.eq."earth") then 175 ierr = NF_DEF_DIM (nid, "index", length, idim_index) 176 ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu) 177 ierr = NF_DEF_DIM (nid, "rlatu", jjp1, idim_rlatu) 178 ierr = NF_DEF_DIM (nid, "rlonv", iip1, idim_rlonv) 179 ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv) 180 ierr = NF_DEF_DIM (nid, "sigs", llm, idim_s) 181 ierr = NF_DEF_DIM (nid, "sig", llmp1, idim_sig) 182 ierr = NF_DEF_DIM (nid, "temps", NF_UNLIMITED, idim_tim) 183 else 184 ierr = NF_DEF_DIM (nid, "index", length, idim_index) 185 ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu) 186 ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu) 187 ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv) 188 ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv) 189 ierr = NF_DEF_DIM (nid, "altitude", llm, idim_s) 190 ierr = NF_DEF_DIM (nid, "interlayer", llmp1, idim_sig) 191 ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_tim) 192 endif 193 c 194 ierr = NF_ENDDEF(nid) ! sortir du mode de definition 195 c 196 c Definir et enregistrer certains champs invariants: 197 c 198 ierr = NF_REDEF (nid) 199 cIM 220306 BEG 200 #ifdef NC_DOUBLE 201 ierr = NF_DEF_VAR (nid,"controle",NF_DOUBLE,1,idim_index,nvarid) 202 #else 203 ierr = NF_DEF_VAR (nid,"controle",NF_FLOAT,1,idim_index,nvarid) 204 #endif 205 cIM 220306 END 206 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, 207 . "Parametres de controle") 208 ierr = NF_ENDDEF(nid) 209 call NF95_PUT_VAR(nid,nvarid,tab_cntrl) 210 c 211 ierr = NF_REDEF (nid) 212 cIM 220306 BEG 213 #ifdef NC_DOUBLE 214 ierr = NF_DEF_VAR (nid,"rlonu",NF_DOUBLE,1,idim_rlonu,nvarid) 215 #else 216 ierr = NF_DEF_VAR (nid,"rlonu",NF_FLOAT,1,idim_rlonu,nvarid) 217 #endif 218 cIM 220306 END 219 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23, 220 . "Longitudes des points U") 221 ierr = NF_ENDDEF(nid) 222 call NF95_PUT_VAR(nid,nvarid,rlonu) 223 c 224 ierr = NF_REDEF (nid) 225 cIM 220306 BEG 226 #ifdef NC_DOUBLE 227 ierr = NF_DEF_VAR (nid,"rlatu",NF_DOUBLE,1,idim_rlatu,nvarid) 228 #else 229 ierr = NF_DEF_VAR (nid,"rlatu",NF_FLOAT,1,idim_rlatu,nvarid) 230 #endif 231 cIM 220306 END 232 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, 233 . "Latitudes des points U") 234 ierr = NF_ENDDEF(nid) 235 call NF95_PUT_VAR (nid,nvarid,rlatu) 236 c 237 ierr = NF_REDEF (nid) 238 cIM 220306 BEG 239 #ifdef NC_DOUBLE 240 ierr = NF_DEF_VAR (nid,"rlonv",NF_DOUBLE,1,idim_rlonv,nvarid) 241 #else 242 ierr = NF_DEF_VAR (nid,"rlonv",NF_FLOAT,1,idim_rlonv,nvarid) 243 #endif 244 cIM 220306 END 245 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23, 246 . "Longitudes des points V") 247 ierr = NF_ENDDEF(nid) 248 call NF95_PUT_VAR(nid,nvarid,rlonv) 249 c 250 ierr = NF_REDEF (nid) 251 cIM 220306 BEG 252 #ifdef NC_DOUBLE 253 ierr = NF_DEF_VAR (nid,"rlatv",NF_DOUBLE,1,idim_rlatv,nvarid) 254 #else 255 ierr = NF_DEF_VAR (nid,"rlatv",NF_FLOAT,1,idim_rlatv,nvarid) 256 #endif 257 cIM 220306 END 258 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, 259 . "Latitudes des points V") 260 ierr = NF_ENDDEF(nid) 261 call NF95_PUT_VAR(nid,nvarid,rlatv) 262 c 263 if (start_file_type.eq."earth") then 264 ierr = NF_REDEF (nid) 265 cIM 220306 BEG 266 #ifdef NC_DOUBLE 267 ierr = NF_DEF_VAR (nid,"nivsigs",NF_DOUBLE,1,idim_s,nvarid) 268 #else 269 ierr = NF_DEF_VAR (nid,"nivsigs",NF_FLOAT,1,idim_s,nvarid) 270 #endif 271 cIM 220306 END 272 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 28, 273 . "Numero naturel des couches s") 274 ierr = NF_ENDDEF(nid) 275 call NF95_PUT_VAR(nid,nvarid,nivsigs) 276 c 277 ierr = NF_REDEF (nid) 278 cIM 220306 BEG 279 #ifdef NC_DOUBLE 280 ierr = NF_DEF_VAR (nid,"nivsig",NF_DOUBLE,1,idim_sig,nvarid) 281 #else 282 ierr = NF_DEF_VAR (nid,"nivsig",NF_FLOAT,1,idim_sig,nvarid) 283 #endif 284 cIM 220306 END 285 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 32, 286 . "Numero naturel des couches sigma") 287 ierr = NF_ENDDEF(nid) 288 call NF95_PUT_VAR(nid,nvarid,nivsig) 289 endif ! of if (start_file_type.eq."earth") 290 c 291 ierr = NF_REDEF (nid) 292 cIM 220306 BEG 293 #ifdef NC_DOUBLE 294 ierr = NF_DEF_VAR (nid,"ap",NF_DOUBLE,1,idim_sig,nvarid) 295 #else 296 ierr = NF_DEF_VAR (nid,"ap",NF_FLOAT,1,idim_sig,nvarid) 297 #endif 298 cIM 220306 END 299 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26, 300 . "Coefficient A pour hybride") 301 ierr = NF_ENDDEF(nid) 302 call NF95_PUT_VAR(nid,nvarid,ap) 303 c 304 ierr = NF_REDEF (nid) 305 cIM 220306 BEG 306 #ifdef NC_DOUBLE 307 ierr = NF_DEF_VAR (nid,"bp",NF_DOUBLE,1,idim_sig,nvarid) 308 #else 309 ierr = NF_DEF_VAR (nid,"bp",NF_FLOAT,1,idim_sig,nvarid) 310 #endif 311 cIM 220306 END 312 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26, 313 . "Coefficient B pour hybride") 314 ierr = NF_ENDDEF(nid) 315 call NF95_PUT_VAR(nid,nvarid,bp) 316 c 317 if (start_file_type.ne."earth") then 318 ierr = NF_REDEF (nid) 319 cIM 220306 BEG 320 #ifdef NC_DOUBLE 321 ierr = NF_DEF_VAR (nid,"aps",NF_DOUBLE,1,idim_s,nvarid) 322 #else 323 ierr = NF_DEF_VAR (nid,"aps",NF_FLOAT,1,idim_s,nvarid) 324 #endif 325 cIM 220306 END 326 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 37, 327 . "Coef AS: hybrid pressure at midlayers") 328 ierr = NF_ENDDEF(nid) 329 call NF95_PUT_VAR(nid,nvarid,aps) 330 c 331 ierr = NF_REDEF (nid) 332 cIM 220306 BEG 333 #ifdef NC_DOUBLE 334 ierr = NF_DEF_VAR (nid,"bps",NF_DOUBLE,1,idim_s,nvarid) 335 #else 336 ierr = NF_DEF_VAR (nid,"bps",NF_FLOAT,1,idim_s,nvarid) 337 #endif 338 cIM 220306 END 339 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 34, 340 . "Coef BS: hybrid sigma at midlayers") 341 ierr = NF_ENDDEF(nid) 342 call NF95_PUT_VAR(nid,nvarid,bps) 343 endif ! of if (start_file_type.ne."earth") 344 c 345 ierr = NF_REDEF (nid) 346 cIM 220306 BEG 347 #ifdef NC_DOUBLE 348 ierr = NF_DEF_VAR (nid,"presnivs",NF_DOUBLE,1,idim_s,nvarid) 349 #else 350 ierr = NF_DEF_VAR (nid,"presnivs",NF_FLOAT,1,idim_s,nvarid) 351 #endif 352 cIM 220306 END 353 ierr = NF_ENDDEF(nid) 354 call NF95_PUT_VAR(nid,nvarid,presnivs) 355 c 356 if (start_file_type.ne."earth") then 76 !!! AS: idecal is a hack to be able to read planeto starts... 77 !!! .... while keeping everything OK for LMDZ EARTH 78 if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then 79 write(lunout,*) trim(modname),' : Planeto-like start file' 80 start_file_type="planeto" 81 idecal = 4 82 else 83 write(lunout,*) trim(modname),' : Earth-like start file' 84 idecal = 5 85 endif 86 87 tab_cntrl(:) = 0. 88 tab_cntrl(1) = REAL(iim) 89 tab_cntrl(2) = REAL(jjm) 90 tab_cntrl(3) = REAL(llm) 91 if (start_file_type.eq."earth") then 92 tab_cntrl(4)=REAL(day_ref) 93 else 94 !tab_cntrl(4)=REAL(day_end) 95 tab_cntrl(4)=REAL(iday_end) 96 endif 97 tab_cntrl(5) = REAL(annee_ref) 98 tab_cntrl(idecal+1) = rad 99 tab_cntrl(idecal+2) = omeg 100 tab_cntrl(idecal+3) = g 101 tab_cntrl(idecal+4) = cpp 102 tab_cntrl(idecal+5) = kappa 103 tab_cntrl(idecal+6) = daysec 104 tab_cntrl(idecal+7) = dtvr 105 tab_cntrl(idecal+8) = etot0 106 tab_cntrl(idecal+9) = ptot0 107 tab_cntrl(idecal+10) = ztot0 108 tab_cntrl(idecal+11) = stot0 109 tab_cntrl(idecal+12) = ang0 110 tab_cntrl(idecal+13) = pa 111 tab_cntrl(idecal+14) = preff 112 113 ! ..... parameters for the zoom ...... 114 tab_cntrl(idecal+15) = clon 115 tab_cntrl(idecal+16) = clat 116 tab_cntrl(idecal+17) = grossismx 117 tab_cntrl(idecal+18) = grossismy 118 ! 119 IF ( fxyhypb ) THEN 120 tab_cntrl(idecal+19) = 1. 121 tab_cntrl(idecal+20) = dzoomx 122 tab_cntrl(idecal+21) = dzoomy 123 tab_cntrl(idecal+22) = 0. 124 tab_cntrl(idecal+23) = taux 125 tab_cntrl(idecal+24) = tauy 126 ELSE 127 tab_cntrl(idecal+19) = 0. 128 tab_cntrl(idecal+20) = dzoomx 129 tab_cntrl(idecal+21) = dzoomy 130 tab_cntrl(idecal+22) = 0. 131 tab_cntrl(idecal+23) = 0. 132 tab_cntrl(idecal+24) = 0. 133 IF( ysinus ) tab_cntrl(idecal+22) = 1. 134 ENDIF 135 136 if (start_file_type.eq."earth") then 137 tab_cntrl(idecal+25) = REAL(iday_end) 138 tab_cntrl(idecal+26) = REAL(itau_dyn + itaufin) 139 ! start_time: start_time of simulation (not necessarily 0.) 140 tab_cntrl(idecal+27) = start_time 141 endif 142 143 if (planet_type=="mars") then ! For Mars only 144 tab_cntrl(29)=hour_ini 145 endif 146 147 !--- File creation 148 CALL err(NF90_CREATE(fichnom,NF90_CLOBBER,nid)) 149 150 !--- Some global attributes 151 CALL err(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier demarrage dynamique")) 152 153 !--- Dimensions 154 if (start_file_type.eq."earth") then 155 CALL err(NF90_DEF_DIM(nid,"index", length, indexID)) 156 CALL err(NF90_DEF_DIM(nid,"rlonu", iip1, rlonuID)) 157 CALL err(NF90_DEF_DIM(nid,"rlatu", jjp1, rlatuID)) 158 CALL err(NF90_DEF_DIM(nid,"rlonv", iip1, rlonvID)) 159 CALL err(NF90_DEF_DIM(nid,"rlatv", jjm, rlatvID)) 160 CALL err(NF90_DEF_DIM(nid,"sigs", llm, sID)) 161 CALL err(NF90_DEF_DIM(nid,"sig", llmp1, sigID)) 162 CALL err(NF90_DEF_DIM(nid,"temps", NF90_UNLIMITED, timID)) 163 else 164 CALL err(NF90_DEF_DIM(nid,"index", length, indexID)) 165 CALL err(NF90_DEF_DIM(nid,"rlonu", iip1, rlonuID)) 166 CALL err(NF90_DEF_DIM(nid,"latitude", jjp1, rlatuID)) 167 CALL err(NF90_DEF_DIM(nid,"longitude", iip1, rlonvID)) 168 CALL err(NF90_DEF_DIM(nid,"rlatv", jjm, rlatvID)) 169 CALL err(NF90_DEF_DIM(nid,"altitude", llm, sID)) 170 CALL err(NF90_DEF_DIM(nid,"interlayer", llmp1, sigID)) 171 CALL err(NF90_DEF_DIM(nid,"Time", NF90_UNLIMITED, timID)) 172 endif 173 174 !--- Define and save invariant fields 175 CALL put_var1(nid,"controle","Parametres de controle" ,[indexID],tab_cntrl) 176 CALL put_var1(nid,"rlonu" ,"Longitudes des points U",[rlonuID],rlonu) 177 CALL put_var1(nid,"rlatu" ,"Latitudes des points U" ,[rlatuID],rlatu) 178 CALL put_var1(nid,"rlonv" ,"Longitudes des points V",[rlonvID],rlonv) 179 CALL put_var1(nid,"rlatv" ,"Latitudes des points V" ,[rlatvID],rlatv) 180 if (start_file_type.eq."earth") then 181 CALL put_var1(nid,"nivsigs" ,"Numero naturel des couches s" ,[sID] ,nivsigs) 182 CALL put_var1(nid,"nivsig" ,"Numero naturel des couches sigma",[sigID],nivsig) 183 endif ! of if (start_file_type.eq."earth") 184 CALL put_var1(nid,"ap" ,"Coefficient A pour hybride" ,[sigID],ap) 185 CALL put_var1(nid,"bp" ,"Coefficient B pour hybride" ,[sigID],bp) 186 if (start_file_type.ne."earth") then 187 CALL put_var1(nid,"aps","Coef AS: hybrid pressure at midlayers",[sID],aps) 188 CALL put_var1(nid,"bps","Coef BS: hybrid sigma at midlayers",[sID],bps) 189 endif ! of if (start_file_type.eq."earth") 190 CALL put_var1(nid,"presnivs","" ,[sID] ,presnivs) 191 if (start_file_type.ne."earth") then 357 192 ierr = NF_REDEF (nid) 358 193 #ifdef NC_DOUBLE 359 ierr = NF_DEF_VAR(nid,"latitude",NF_DOUBLE,1, idim_rlatu,nvarid)194 ierr = NF_DEF_VAR(nid,"latitude",NF_DOUBLE,1,rlatuID,vID) 360 195 #else 361 ierr = NF_DEF_VAR(nid,"latitude",NF_FLOAT,1, idim_rlatu,nvarid)196 ierr = NF_DEF_VAR(nid,"latitude",NF_FLOAT,1,rlatuID,vID) 362 197 #endif 363 ierr =NF_PUT_ATT_TEXT(nid, nvarid,'units',13,"degrees_north")364 ierr = NF_PUT_ATT_TEXT (nid, nvarid,"long_name", 14,365 ."North latitude")198 ierr =NF_PUT_ATT_TEXT(nid,vID,'units',13,"degrees_north") 199 ierr = NF_PUT_ATT_TEXT (nid,vID,"long_name", 14, & 200 "North latitude") 366 201 ierr = NF_ENDDEF(nid) 367 call NF95_PUT_VAR(nid, nvarid,rlatu*180/pi)368 c 202 call NF95_PUT_VAR(nid,vID,rlatu*180/pi) 203 ! 369 204 ierr = NF_REDEF (nid) 370 205 #ifdef NC_DOUBLE 371 ierr=NF_DEF_VAR(nid,"longitude",NF_DOUBLE,1, idim_rlonv,nvarid)206 ierr=NF_DEF_VAR(nid,"longitude",NF_DOUBLE,1,rlonvID,vID) 372 207 #else 373 ierr=NF_DEF_VAR(nid,"longitude",NF_FLOAT,1, idim_rlonv,nvarid)208 ierr=NF_DEF_VAR(nid,"longitude",NF_FLOAT,1,rlonvID,vID) 374 209 #endif 375 ierr = NF_PUT_ATT_TEXT (nid, nvarid,"long_name", 14,376 ."East longitude")377 ierr = NF_PUT_ATT_TEXT(nid, nvarid,'units',12,"degrees_east")210 ierr = NF_PUT_ATT_TEXT (nid,vID,"long_name", 14, & 211 "East longitude") 212 ierr = NF_PUT_ATT_TEXT(nid,vID,'units',12,"degrees_east") 378 213 ierr = NF_ENDDEF(nid) 379 call NF95_PUT_VAR(nid, nvarid,rlonv*180/pi)380 c 214 call NF95_PUT_VAR(nid,vID,rlonv*180/pi) 215 ! 381 216 ierr = NF_REDEF (nid) 382 217 #ifdef NC_DOUBLE 383 ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, 384 . idim_s,nvarid)218 ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, & 219 sID,vID) 385 220 #else 386 ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, 387 . idim_s,nvarid)221 ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, & 222 sID,vID) 388 223 #endif 389 ierr = NF_PUT_ATT_TEXT(nid, nvarid,"long_name",10,"pseudo-alt")390 ierr = NF_PUT_ATT_TEXT (nid, nvarid,'units',2,"km")391 ierr = NF_PUT_ATT_TEXT (nid, nvarid,'positive',2,"up")224 ierr = NF_PUT_ATT_TEXT(nid,vID,"long_name",10,"pseudo-alt") 225 ierr = NF_PUT_ATT_TEXT (nid,vID,'units',2,"km") 226 ierr = NF_PUT_ATT_TEXT (nid,vID,'positive',2,"up") 392 227 ierr = NF_ENDDEF(nid) 393 call NF95_PUT_VAR(nid,nvarid,pseudoalt) 394 endif ! of if (start_file_type.ne."earth") 395 c 396 c Coefficients de passage cov. <-> contra. <--> naturel 397 c 398 ierr = NF_REDEF (nid) 399 dims2(1) = idim_rlonu 400 dims2(2) = idim_rlatu 401 cIM 220306 BEG 402 #ifdef NC_DOUBLE 403 ierr = NF_DEF_VAR (nid,"cu",NF_DOUBLE,2,dims2,nvarid) 404 #else 405 ierr = NF_DEF_VAR (nid,"cu",NF_FLOAT,2,dims2,nvarid) 406 #endif 407 cIM 220306 END 408 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29, 409 . "Coefficient de passage pour U") 410 ierr = NF_ENDDEF(nid) 411 call NF95_PUT_VAR(nid,nvarid,cu) 412 c 413 ierr = NF_REDEF (nid) 414 dims2(1) = idim_rlonv 415 dims2(2) = idim_rlatv 416 cIM 220306 BEG 417 #ifdef NC_DOUBLE 418 ierr = NF_DEF_VAR (nid,"cv",NF_DOUBLE,2,dims2,nvarid) 419 #else 420 ierr = NF_DEF_VAR (nid,"cv",NF_FLOAT,2,dims2,nvarid) 421 #endif 422 cIM 220306 END 423 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29, 424 . "Coefficient de passage pour V") 425 ierr = NF_ENDDEF(nid) 426 call NF95_PUT_VAR(nid,nvarid,cv) 427 c 428 c Aire de chaque maille: 429 c 430 ierr = NF_REDEF (nid) 431 dims2(1) = idim_rlonv 432 dims2(2) = idim_rlatu 433 cIM 220306 BEG 434 #ifdef NC_DOUBLE 435 ierr = NF_DEF_VAR (nid,"aire",NF_DOUBLE,2,dims2,nvarid) 436 #else 437 ierr = NF_DEF_VAR (nid,"aire",NF_FLOAT,2,dims2,nvarid) 438 #endif 439 cIM 220306 END 440 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22, 441 . "Aires de chaque maille") 442 ierr = NF_ENDDEF(nid) 443 call NF95_PUT_VAR(nid,nvarid,aire) 444 c 445 c Geopentiel au sol: 446 c 447 ierr = NF_REDEF (nid) 448 dims2(1) = idim_rlonv 449 dims2(2) = idim_rlatu 450 cIM 220306 BEG 451 #ifdef NC_DOUBLE 452 ierr = NF_DEF_VAR (nid,"phisinit",NF_DOUBLE,2,dims2,nvarid) 453 #else 454 ierr = NF_DEF_VAR (nid,"phisinit",NF_FLOAT,2,dims2,nvarid) 455 #endif 456 cIM 220306 END 457 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19, 458 . "Geopotentiel au sol") 459 ierr = NF_ENDDEF(nid) 460 call NF95_PUT_VAR(nid,nvarid,phis) 461 c 462 c Definir les variables pour pouvoir les enregistrer plus tard: 463 c 464 ierr = NF_REDEF (nid) ! entrer dans le mode de definition 465 c 466 if (start_file_type.eq."earth") then 467 cIM 220306 BEG 468 #ifdef NC_DOUBLE 469 ierr = NF_DEF_VAR (nid,"temps",NF_DOUBLE,1,idim_tim,nvarid) 470 #else 471 ierr = NF_DEF_VAR (nid,"temps",NF_FLOAT,1,idim_tim,nvarid) 472 #endif 473 cIM 220306 END 474 else ! start_file_type=="planeto" 475 #ifdef NC_DOUBLE 476 ierr = NF_DEF_VAR (nid,"Time",NF_DOUBLE,1,idim_tim,nvarid) 477 #else 478 ierr = NF_DEF_VAR (nid,"Time",NF_FLOAT,1,idim_tim,nvarid) 479 #endif 480 endif ! of if (start_file_type.eq."earth") 481 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19, 482 . "Temps de simulation") 483 write(unites,200)yyears0,mmois0,jjour0 484 200 format('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00') 485 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "units", 30, 486 . unites) 487 488 c 489 dims4(1) = idim_rlonu 490 dims4(2) = idim_rlatu 491 dims4(3) = idim_s 492 dims4(4) = idim_tim 493 cIM 220306 BEG 494 #ifdef NC_DOUBLE 495 ierr = NF_DEF_VAR (nid,"ucov",NF_DOUBLE,4,dims4,nvarid) 496 #else 497 ierr = NF_DEF_VAR (nid,"ucov",NF_FLOAT,4,dims4,nvarid) 498 #endif 499 cIM 220306 END 500 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9, 501 . "Vitesse U") 502 c 503 dims4(1) = idim_rlonv 504 dims4(2) = idim_rlatv 505 dims4(3) = idim_s 506 dims4(4) = idim_tim 507 cIM 220306 BEG 508 #ifdef NC_DOUBLE 509 ierr = NF_DEF_VAR (nid,"vcov",NF_DOUBLE,4,dims4,nvarid) 510 #else 511 ierr = NF_DEF_VAR (nid,"vcov",NF_FLOAT,4,dims4,nvarid) 512 #endif 513 cIM 220306 END 514 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9, 515 . "Vitesse V") 516 c 517 dims4(1) = idim_rlonv 518 dims4(2) = idim_rlatu 519 dims4(3) = idim_s 520 dims4(4) = idim_tim 521 cIM 220306 BEG 522 #ifdef NC_DOUBLE 523 ierr = NF_DEF_VAR (nid,"teta",NF_DOUBLE,4,dims4,nvarid) 524 #else 525 ierr = NF_DEF_VAR (nid,"teta",NF_FLOAT,4,dims4,nvarid) 526 #endif 527 cIM 220306 END 528 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 11, 529 . "Temperature") 530 c 531 dims4(1) = idim_rlonv 532 dims4(2) = idim_rlatu 533 dims4(3) = idim_s 534 dims4(4) = idim_tim 535 536 DO iq=1,nqtot 537 cIM 220306 BEG 538 #ifdef NC_DOUBLE 539 ierr = NF_DEF_VAR (nid,tname(iq),NF_DOUBLE,4,dims4,nvarid) 540 #else 541 ierr = NF_DEF_VAR (nid,tname(iq),NF_FLOAT,4,dims4,nvarid) 542 #endif 543 cIM 220306 END 544 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,ttext(iq)) 545 ENDDO 546 c 547 dims4(1) = idim_rlonv 548 dims4(2) = idim_rlatu 549 dims4(3) = idim_s 550 dims4(4) = idim_tim 551 cIM 220306 BEG 552 #ifdef NC_DOUBLE 553 ierr = NF_DEF_VAR (nid,"masse",NF_DOUBLE,4,dims4,nvarid) 554 #else 555 ierr = NF_DEF_VAR (nid,"masse",NF_FLOAT,4,dims4,nvarid) 556 #endif 557 cIM 220306 END 558 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12, 559 . "C est quoi ?") 560 c 561 dims3(1) = idim_rlonv 562 dims3(2) = idim_rlatu 563 dims3(3) = idim_tim 564 cIM 220306 BEG 565 #ifdef NC_DOUBLE 566 ierr = NF_DEF_VAR (nid,"ps",NF_DOUBLE,3,dims3,nvarid) 567 #else 568 ierr = NF_DEF_VAR (nid,"ps",NF_FLOAT,3,dims3,nvarid) 569 #endif 570 cIM 220306 END 571 ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 15, 572 . "Pression au sol") 573 c 574 ierr = NF_ENDDEF(nid) ! sortir du mode de definition 575 ierr = NF_CLOSE(nid) ! fermer le fichier 576 577 PRINT*,'iim,jjm,llm,iday_end',iim,jjm,llm,iday_end 578 PRINT*,'rad,omeg,g,cpp,kappa', 579 , rad,omeg,g,cpp,kappa 580 581 endif ! mpi_rank==0 582 RETURN 583 END 228 call NF95_PUT_VAR(nid,vID,pseudoalt) 229 CALL err(NF_REDEF(nid)) 230 endif ! of if (start_file_type.ne."earth") 231 232 ! covariant <-> contravariant <-> natural conversion coefficients 233 CALL put_var2(nid,"cu","Coefficient de passage pour U",[rlonuID,rlatuID],cu) 234 CALL put_var2(nid,"cv","Coefficient de passage pour V",[rlonvID,rlatvID],cv) 235 CALL put_var2(nid,"aire","Aires de chaque maille" ,[rlonvID,rlatuID],aire) 236 CALL put_var2(nid,"phisinit","Geopotentiel au sol" ,[rlonvID,rlatuID],phis) 237 238 239 ! Define variables that will be stored later: 240 WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')"),& 241 yyears0,mmois0,jjour0 242 IF (planet_type.eq."earth") THEN 243 CALL cre_var(nid,"temps","Temps de simulation",[timID],unites) 244 ELSE 245 CALL cre_var(nid,"Time","Temps de simulation",[timID],unites) 246 ENDIF 247 248 CALL cre_var(nid,"ucov" ,"Vitesse U" ,[rlonuID,rlatuID,sID,timID]) 249 CALL cre_var(nid,"vcov" ,"Vitesse V" ,[rlonvID,rlatvID,sID,timID]) 250 CALL cre_var(nid,"teta" ,"Temperature",[rlonvID,rlatuID,sID,timID]) 251 252 IF(nqtot.GE.1) THEN 253 DO iq=1,nqtot 254 CALL cre_var(nid,tname(iq),ttext(iq),[rlonvID,rlatuID,sID,timID]) 255 END DO 256 ENDIF 257 258 CALL cre_var(nid,"masse","Masse d air" ,[rlonvID,rlatuID,sID,timID]) 259 CALL cre_var(nid,"ps" ,"Pression au sol",[rlonvID,rlatuID ,timID]) 260 261 CALL err(NF90_CLOSE (nid)) ! close file 262 263 WRITE(lunout,*)TRIM(modname)//': iim,jjm,llm,iday_end',iim,jjm,llm,iday_end 264 WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa 265 266 endif ! of if (mpi_rank==0) 267 268 END SUBROUTINE dynredem0_p 584 269 585 270 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 586 271 587 SUBROUTINE dynredem1_p(fichnom,time, 588 . vcov,ucov,teta,q,masse,ps) 589 USE parallel_lmdz 590 USE infotrac 591 USE control_mod, only : planet_type 592 use netcdf, only: NF90_get_VAR 593 use netcdf95, only: NF95_PUT_VAR 594 USE temps_mod, ONLY: itau_dyn,itaufin 595 596 IMPLICIT NONE 597 c================================================================= 598 c Ecriture du fichier de redemarrage sous format NetCDF 599 c================================================================= 600 #include "dimensions.h" 601 #include "paramet.h" 602 #include "netcdf.inc" 603 #include "comgeom.h" 604 #include "iniprint.h" 605 606 INTEGER l 607 REAL vcov(iip1,jjm,llm),ucov(iip1, jjp1,llm) 608 REAL teta(iip1, jjp1,llm) 609 REAL ps(iip1, jjp1),masse(iip1, jjp1,llm) 610 REAL q(iip1, jjp1, llm, nqtot) 611 CHARACTER*(*) fichnom 612 613 REAL time 614 INTEGER nid, nvarid, nid_trac, nvarid_trac 615 REAL trac_tmp(ip1jmp1,llm) 616 INTEGER ierr, ierr_file 617 INTEGER iq 618 INTEGER length 619 PARAMETER (length = 100) 620 REAL tab_cntrl(length) ! tableau des parametres du run 621 character(len=*),parameter :: modname='dynredem1' 622 character*80 abort_message 623 c 624 INTEGER nb 625 SAVE nb 626 DATA nb / 0 / 627 628 logical exist_file 629 character(len=12) :: start_file_type="earth" ! default start file type 630 631 if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then 632 write(lunout,*) trim(modname),' : Planeto-like start file' 633 start_file_type="planeto" 634 else 635 write(lunout,*) trim(modname),' : Earth-like start file' 636 endif 637 638 call Gather_Field(ucov,ip1jmp1,llm,0) 639 call Gather_Field(vcov,ip1jm,llm,0) 640 call Gather_Field(teta,ip1jmp1,llm,0) 641 call Gather_Field(masse,ip1jmp1,llm,0) 642 call Gather_Field(ps,ip1jmp1,1,0) 272 SUBROUTINE dynredem1_p(fichnom,time,vcov,ucov,teta,q,masse,ps) 273 ! 274 !------------------------------------------------------------------------------- 275 ! Purpose: Write the NetCDF restart file (append). 276 !------------------------------------------------------------------------------- 277 USE parallel_lmdz, ONLY: mpi_rank, gather_field 278 USE infotrac, ONLY: nqtot, tname, type_trac 279 USE control_mod, only : planet_type 280 USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_GET_VAR, NF90_INQ_VARID, & 281 NF90_CLOSE, NF90_WRITE, NF90_PUT_VAR, NF90_NoErr 282 use netcdf95, only: NF95_PUT_VAR 283 USE temps_mod, ONLY: itaufin,itau_dyn 284 USE dynredem_mod, ONLY: dynredem_write_u, dynredem_write_v, dynredem_read_u, & 285 err, modname, fil, msg 286 287 IMPLICIT NONE 288 include "dimensions.h" 289 include "paramet.h" 290 include "netcdf.inc" 291 include "comgeom.h" 292 include "iniprint.h" 293 !=============================================================================== 294 ! Arguments: 295 CHARACTER(LEN=*), INTENT(IN) :: fichnom !-- FILE NAME 296 REAL, INTENT(IN) :: time !-- TIME 297 REAL, INTENT(IN) :: vcov(iip1,jjm, llm) !-- V COVARIANT WIND 298 REAL, INTENT(IN) :: ucov(iip1,jjp1,llm) !-- U COVARIANT WIND 299 REAL, INTENT(IN) :: teta(iip1,jjp1,llm) !-- POTENTIAL TEMPERATURE 300 REAL, INTENT(INOUT) :: q(iip1,jjp1,llm,nqtot) !-- TRACERS 301 REAL, INTENT(IN) :: masse(iip1,jjp1,llm) !-- MASS PER CELL 302 REAL, INTENT(IN) :: ps(iip1,jjp1) !-- GROUND PRESSURE 303 !=============================================================================== 304 ! Local variables: 305 INTEGER :: l, iq, nid, vID, ierr, nid_trac, vID_trac 306 INTEGER,SAVE :: nb=0 307 INTEGER, PARAMETER :: length=100 308 REAL :: tab_cntrl(length) ! tableau des parametres du run 309 CHARACTER(LEN=256) :: var, dum 310 LOGICAL :: lread_inca 311 CHARACTER(LEN=80) :: abort_message 312 CHARACTER(len=12) :: start_file_type="earth" ! default start file type 313 314 ! fill dynredem_mod module variables 315 modname='dynredem1_p'; fil=fichnom 316 317 ! Gather datasets 318 call Gather_Field(ucov,ip1jmp1,llm,0) 319 call Gather_Field(vcov,ip1jm,llm,0) 320 call Gather_Field(teta,ip1jmp1,llm,0) 321 call Gather_Field(masse,ip1jmp1,llm,0) 322 call Gather_Field(ps,ip1jmp1,1,0) 643 323 644 do iq=1,nqtot 645 call Gather_Field(q(:,:,:,iq),ip1jmp1,llm,0) 646 enddo 647 648 649 if (mpi_rank==0) then 650 651 ierr = NF_OPEN(fichnom, NF_WRITE, nid) 652 IF (ierr .NE. NF_NOERR) THEN 653 PRINT*, "dynredem1: Pb. d ouverture "//trim(fichnom) 654 CALL abort 655 ENDIF 656 657 c Ecriture/extension de la coordonnee temps 658 659 nb = nb + 1 660 if (start_file_type.eq."earth") then 661 ierr = NF_INQ_VARID(nid, "temps", nvarid) 324 do iq=1,nqtot 325 call Gather_Field(q(:,:,:,iq),ip1jmp1,llm,0) 326 enddo 327 328 IF (mpi_rank==0) THEN ! only the master writes restart file 329 330 if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then 331 write(lunout,*) trim(modname),' : Planeto-like start file' 332 start_file_type="planeto" 333 else 334 write(lunout,*) trim(modname),' : Earth-like start file' 335 endif 336 337 CALL err(NF90_OPEN(fil,NF90_WRITE,nid),"open",fil) 338 339 !--- Write/extend time coordinate 340 nb = nb + 1 341 if (start_file_type.eq."earth") then 342 ierr = NF_INQ_VARID(nid, "temps", vID) 662 343 IF (ierr .NE. NF_NOERR) THEN 663 344 write(lunout,*) NF_STRERROR(ierr) … … 665 346 CALL abort_gcm(modname,abort_message,ierr) 666 347 ENDIF 667 668 ierr = NF_INQ_VARID(nid,"Time", nvarid)348 else 349 ierr = NF_INQ_VARID(nid,"Time", vID) 669 350 IF (ierr .NE. NF_NOERR) THEN 670 351 write(lunout,*) NF_STRERROR(ierr) … … 672 353 CALL abort_gcm(modname,abort_message,ierr) 673 354 ENDIF 674 endif ! of if (start_file_type.eq."earth") 675 call NF95_PUT_VAR(nid,nvarid,time,start=(/nb/)) 676 PRINT*, "Enregistrement pour ", nb, time 677 678 c 679 c Re-ecriture du tableau de controle, itaufin n'est plus defini quand 680 c on passe dans dynredem0 681 ierr = NF_INQ_VARID (nid, "controle", nvarid) 682 IF (ierr .NE. NF_NOERR) THEN 683 abort_message="dynredem1: Le champ <controle> est absent" 684 ierr = 1 685 CALL abort_gcm(modname,abort_message,ierr) 686 ENDIF 687 ierr = NF90_GET_VAR(nid, nvarid, tab_cntrl) 688 if (start_file_type=="earth") then 689 tab_cntrl(31) = REAL(itau_dyn + itaufin) 690 else 691 tab_cntrl(31) = 0 692 endif 693 call NF95_PUT_VAR(nid,nvarid,tab_cntrl) 694 695 c Ecriture des champs 696 c 697 ierr = NF_INQ_VARID(nid, "ucov", nvarid) 698 IF (ierr .NE. NF_NOERR) THEN 699 PRINT*, "Variable ucov n est pas definie" 700 CALL abort 701 ENDIF 702 call NF95_PUT_VAR(nid,nvarid,ucov,start=(/1,1,1,nb/)) 703 704 ierr = NF_INQ_VARID(nid, "vcov", nvarid) 705 IF (ierr .NE. NF_NOERR) THEN 706 PRINT*, "Variable vcov n est pas definie" 707 CALL abort 708 ENDIF 709 call NF95_PUT_VAR(nid,nvarid,vcov,start=(/1,1,1,nb/)) 710 711 ierr = NF_INQ_VARID(nid, "teta", nvarid) 712 IF (ierr .NE. NF_NOERR) THEN 713 PRINT*, "Variable teta n est pas definie" 714 CALL abort 715 ENDIF 716 call NF95_PUT_VAR(nid,nvarid,teta,start=(/1,1,1,nb/)) 717 718 IF (type_trac == 'inca') THEN 719 ! Ajout Anne pour lecture valeurs traceurs dans un fichier start_trac.nc 720 inquire(FILE="start_trac.nc", EXIST=exist_file) 721 print *, "EXIST", exist_file 722 if (exist_file) then 723 ierr_file = NF_OPEN ("start_trac.nc", NF_NOWRITE,nid_trac) 724 IF (ierr_file .NE.NF_NOERR) THEN 725 write(6,*)' Pb d''ouverture du fichier start_trac.nc' 726 write(6,*)' ierr = ', ierr_file 727 ENDIF 728 else 729 ierr_file = 2 730 endif 731 END IF 732 733 do iq=1,nqtot 734 735 IF (type_trac /= 'inca') THEN 736 ierr = NF_INQ_VARID(nid, tname(iq), nvarid) 737 IF (ierr .NE. NF_NOERR) THEN 738 PRINT*, "Variable tname(iq) n est pas definie" 739 CALL abort 740 ENDIF 741 call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq),start=(/1,1,1,nb/)) 742 ELSE ! type_trac = inca 743 ! lecture de la valeur du traceur dans start_trac.nc 744 IF (ierr_file .ne. 2) THEN 745 ierr = NF_INQ_VARID (nid_trac, tname(iq), nvarid_trac) 746 IF (ierr .NE. NF_NOERR) THEN 747 PRINT*, tname(iq),"est absent de start_trac.nc" 748 ierr = NF_INQ_VARID(nid, tname(iq), nvarid) 749 IF (ierr .NE. NF_NOERR) THEN 750 PRINT*, "Variable ", tname(iq)," n est pas definie" 751 CALL abort 752 ENDIF 753 call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq)) 754 755 ELSE 756 PRINT*, tname(iq), "est present dans start_trac.nc" 757 ierr = NF90_GET_VAR(nid_trac, nvarid_trac, trac_tmp) 758 IF (ierr .NE. NF_NOERR) THEN 759 PRINT*, "Lecture echouee pour", tname(iq) 760 CALL abort 761 ENDIF 762 ierr = NF_INQ_VARID(nid, tname(iq), nvarid) 763 IF (ierr .NE. NF_NOERR) THEN 764 PRINT*, "Variable ", tname(iq)," n est pas definie" 765 CALL abort 766 ENDIF 767 call NF95_PUT_VAR(nid, nvarid, trac_tmp) 768 769 ENDIF ! IF (ierr .NE. NF_NOERR) 770 ! fin lecture du traceur 771 ELSE ! si il n'y a pas de fichier start_trac.nc 772 ! print *, 'il n y a pas de fichier start_trac' 773 ierr = NF_INQ_VARID(nid, tname(iq), nvarid) 774 IF (ierr .NE. NF_NOERR) THEN 775 PRINT*, "Variable tname(iq) n est pas definie" 776 CALL abort 777 ENDIF 778 call NF95_PUT_VAR(nid,nvarid,q(:,:,:,iq), 779 & start=(/1,1,1,nb/)) 780 ENDIF ! (ierr_file .ne. 2) 781 END IF !type_trac 782 783 ENDDO 784 785 786 787 c 788 ierr = NF_INQ_VARID(nid, "masse", nvarid) 789 IF (ierr .NE. NF_NOERR) THEN 790 PRINT*, "Variable masse n est pas definie" 791 CALL abort 792 ENDIF 793 call NF95_PUT_VAR(nid,nvarid,masse,start=(/1,1,1,nb/)) 794 c 795 ierr = NF_INQ_VARID(nid, "ps", nvarid) 796 IF (ierr .NE. NF_NOERR) THEN 797 PRINT*, "Variable ps n est pas definie" 798 CALL abort 799 ENDIF 800 call NF95_PUT_VAR(nid,nvarid,ps,start=(/1,1,nb/)) 801 802 ierr = NF_CLOSE(nid) 803 c 804 endif ! mpi_rank==0 805 806 RETURN 807 END 808 355 endif ! of if (start_file_type.eq."earth") 356 call NF95_PUT_VAR(nid,vID,time,start=(/nb/)) 357 WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time 358 359 !--- Rewrite control table (itaufin undefined in dynredem0) 360 var="controle" 361 CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var) 362 CALL err(NF90_GET_VAR(nid,vID,tab_cntrl),"get",var) 363 if (start_file_type=="earth") then 364 tab_cntrl(31) = REAL(itau_dyn + itaufin) 365 else 366 tab_cntrl(31) = 0 367 endif 368 CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var) 369 CALL err(NF90_PUT_VAR(nid,vID,tab_cntrl),"put",var) 370 371 !--- Save fields 372 CALL dynredem_write_u(nid,"ucov" ,ucov ,llm, nb) 373 CALL dynredem_write_v(nid,"vcov" ,vcov ,llm, nb) 374 CALL dynredem_write_u(nid,"teta" ,teta ,llm, nb) 375 CALL dynredem_write_u(nid,"masse",masse,llm, nb) 376 CALL dynredem_write_u(nid,"ps" ,ps ,1, nb) 377 378 !--- Tracers in file "start_trac.nc" (added by Anne) 379 lread_inca=.FALSE.; fil="start_trac.nc" 380 IF(type_trac=='inca') INQUIRE(FILE=fil,EXIST=lread_inca) 381 IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open") 382 383 !--- Save tracers 384 IF(nqtot.GE.1) THEN 385 DO iq=1,nqtot 386 var=tname(iq); ierr=-1 387 IF(lread_inca) THEN !--- Possibly read from "start_trac.nc" 388 fil="start_trac.nc" 389 ierr=NF90_INQ_VARID(nid_trac,var,vID_trac) 390 dum='inq'; IF(ierr==NF90_NoErr) dum='fnd' 391 WRITE(lunout,*)msg(dum,var) 392 393 394 IF(ierr==NF90_NoErr) CALL dynredem_read_u(nid_trac,var,q(:,:,:,iq),llm) 395 END IF ! of IF(lread_inca) 396 fil=fichnom 397 CALL dynredem_write_u(nid,var,q(:,:,:,iq),llm,nb) 398 END DO ! of DO iq=1,nqtot 399 ENDIF ! of IF(nqtot.GE.1) 400 401 CALL err(NF90_CLOSE(nid),"close") 402 fil="start_trac.nc" 403 IF(lread_inca) CALL err(NF90_CLOSE(nid_trac),"close") 404 405 ENDIF ! of IF (mpi_rank==0) 406 407 END SUBROUTINE dynredem1_p 408 -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F ¶
r1422 r1508 5 5 $ ( nq,vcovm1,ucovm1,tetam1,psm1,massem1, 6 6 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis) !,finvmaold) 7 USE parallel_lmdz 7 USE parallel_lmdz, ONLY: ij_begin, ij_end, pole_nord, pole_sud 8 8 USE control_mod, only : planet_type 9 9 USE comvert_mod, ONLY: ap,bp -
TabularUnified trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F ¶
r1422 r1508 771 771 ! CALL FTRACE_REGION_BEGIN("integrd") 772 772 773 CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,773 CALL integrd_p (nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 , 774 774 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ) 775 775 ! $ finvmaold ) … … 886 886 IF (ip_ebil_dyn.ge.1 ) THEN 887 887 ztit='bil dyn' 888 ! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!888 ! Ehouarn: be careful, diagedyn is Earth-specific! 889 889 IF (planet_type.eq."earth") THEN 890 #ifdef CPP_EARTH891 890 CALL diagedyn(ztit,2,1,1,dtphys 892 891 & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2)) 893 #endif894 892 ENDIF 895 893 ENDIF -
TabularUnified trunk/LMDZ.COMMON/libf/misc/wxios.F90 ¶
r1441 r1508 21 21 CHARACTER(len=100) :: g_field_name = "nofield" 22 22 !$OMP THREADPRIVATE(g_flag_xml,g_field_name) 23 23 REAL :: missing_val_omp 24 REAL :: missing_val 25 !$OMP THREADPRIVATE(missing_val) 24 26 25 27 CONTAINS
Note: See TracChangeset
for help on using the changeset viewer.
