Changeset 1665 for LMDZ5/branches/testing/libf/dyn3d
- Timestamp:
- Oct 9, 2012, 3:35:26 PM (12 years ago)
- Location:
- LMDZ5/branches/testing
- Files:
-
- 1 deleted
- 21 edited
- 1 copied
-
. (modified) (1 prop)
-
libf/dyn3d/calfis.F (modified) (6 diffs)
-
libf/dyn3d/ce0l.F90 (modified) (1 diff)
-
libf/dyn3d/comvert.h (modified) (2 diffs)
-
libf/dyn3d/conf_gcm.F (modified) (1 diff)
-
libf/dyn3d/control_mod.F90 (modified) (1 diff)
-
libf/dyn3d/disvert.F90 (modified) (4 diffs)
-
libf/dyn3d/dynetat0.F (modified) (1 diff)
-
libf/dyn3d/dynredem.F (modified) (1 diff)
-
libf/dyn3d/etat0_netcdf.F90 (modified) (1 diff)
-
libf/dyn3d/exner_hyb.F (modified) (1 diff)
-
libf/dyn3d/exner_milieu.F (modified) (1 diff)
-
libf/dyn3d/gcm.F (modified) (7 diffs)
-
libf/dyn3d/guide_mod.F90 (modified) (1 diff)
-
libf/dyn3d/iniacademic.F90 (modified) (1 diff)
-
libf/dyn3d/iniconst.F (deleted)
-
libf/dyn3d/iniconst.F90 (copied) (copied from LMDZ5/trunk/libf/dyn3d/iniconst.F90)
-
libf/dyn3d/inidissip.F90 (modified) (3 diffs)
-
libf/dyn3d/inigrads.F (modified) (3 diffs)
-
libf/dyn3d/integrd.F (modified) (7 diffs)
-
libf/dyn3d/leapfrog.F (modified) (17 diffs)
-
libf/dyn3d/temps.h (modified) (1 diff)
-
libf/dyn3d/wrgrads.F (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/branches/testing
- Property svn:mergeinfo changed
/LMDZ5/trunk merged: 1576-1580,1582,1584,1591-1593,1597-1598,1600,1604-1620,1622-1628
- Property svn:mergeinfo changed
-
LMDZ5/branches/testing/libf/dyn3d/calfis.F
r1407 r1665 165 165 PARAMETER(ntetaSTD=3) 166 166 REAL rtetaSTD(ntetaSTD) 167 DATA rtetaSTD/350., 380., 405./ 167 DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !! 168 168 REAL PVteta(ngridmx,ntetaSTD) 169 169 c … … 434 434 c 435 435 if (planet_type=="earth") then 436 #ifdef CPP_EARTH 436 #ifdef CPP_PHYS 437 ! PVtheta calls tetalevel, which is in the physics 437 438 cIM calcul PV a teta=350, 380, 405K 438 439 CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta, … … 450 451 451 452 452 if (planet_type=="earth") then453 #ifdef CPP_EARTH454 453 455 454 ! write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys … … 460 459 zdqfic(:,:,:)=0. 461 460 462 do isplit=1,nsplit_phys 461 if (planet_type=="earth") then 462 #ifdef CPP_PHYS 463 464 do isplit=1,nsplit_phys 463 465 464 466 jH_cur_split=jH_cur+(isplit-1) * dtvr / (daysec *nsplit_phys) … … 503 505 zdqfic(:,:,:)=zdqfic(:,:,:)+zdqfi(:,:,:) 504 506 505 enddo 507 enddo ! of do isplit=1,nsplit_phys 508 509 #endif 510 ! of #ifdef CPP_PHYS 511 endif ! of if (planet_type=="earth") 512 506 513 zdufi(:,:)=zdufic(:,:)/nsplit_phys 507 514 zdvfi(:,:)=zdvfic(:,:)/nsplit_phys … … 509 516 zdqfi(:,:,:)=zdqfic(:,:,:)/nsplit_phys 510 517 511 #endif512 endif !of if (planet_type=="earth")513 518 514 519 500 CONTINUE -
LMDZ5/branches/testing/libf/dyn3d/ce0l.F90
r1664 r1665 28 28 IMPLICIT NONE 29 29 #ifndef CPP_EARTH 30 #include "iniprint.h" 30 31 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' 31 32 #else -
LMDZ5/branches/testing/libf/dyn3d/comvert.h
r1520 r1665 9 9 & aps(llm),bps(llm),scaleheight 10 10 11 common/comverti/disvert_type 11 common/comverti/disvert_type, pressure_exner 12 12 13 13 real ap ! hybrid pressure contribution at interlayers … … 30 30 ! using 'z2sig.def' (or 'esasig.def) file 31 31 32 logical pressure_exner 33 ! compute pressure inside layers using Exner function, else use mean 34 ! of pressure values at interfaces 35 32 36 !----------------------------------------------------------------------- -
LMDZ5/branches/testing/libf/dyn3d/conf_gcm.F
r1664 r1665 155 155 nday = 10 156 156 CALL getin('nday',nday) 157 158 !Config Key = starttime 159 !Config Desc = Heure de depart de la simulation 160 !Config Def = 0 161 !Config Help = Heure de depart de la simulation 162 !Config en jour 163 starttime = 0 164 CALL getin('starttime',starttime) 157 165 158 166 !Config Key = day_step -
LMDZ5/branches/testing/libf/dyn3d/control_mod.F90
r1502 r1665 10 10 IMPLICIT NONE 11 11 12 REAL :: periodav 12 REAL :: periodav, starttime 13 13 INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys 14 14 INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy -
LMDZ5/branches/testing/libf/dyn3d/disvert.F90
r1524 r1665 4 4 5 5 ! Auteur : P. Le Van 6 7 use new_unit_m, only: new_unit 8 use ioipsl, only: getin 6 9 7 10 IMPLICIT NONE … … 26 29 real zk, zkm1, dzk1, dzk2, k0, k1 27 30 28 INTEGER l 31 INTEGER l, unit 29 32 REAL dsigmin 30 33 REAL alpha, beta, deltaz 31 INTEGER iostat32 34 REAL x 33 35 character(len=*),parameter :: modname="disvert" 34 36 37 character(len=6):: vert_sampling 38 ! (allowed values are "param", "tropo", "strato" and "read") 39 35 40 !----------------------------------------------------------------------- 41 42 print *, "Call sequence information: disvert" 36 43 37 44 ! default scaleheight is 8km for earth 38 45 scaleheight=8. 39 46 40 OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat) 47 vert_sampling = merge("strato", "tropo ", ok_strato) ! default value 48 call getin('vert_sampling', vert_sampling) 49 print *, 'vert_sampling = ' // vert_sampling 41 50 42 IF (iostat == 0) THEN 43 ! cas 1 on lit les options dans sigma.def: 51 select case (vert_sampling) 52 case ("param") 53 ! On lit les options dans sigma.def: 54 OPEN(99, file='sigma.def', status='old', form='formatted') 44 55 READ(99, *) scaleheight ! hauteur d'echelle 8. 45 56 READ(99, *) deltaz ! epaiseur de la premiere couche 0.04 … … 68 79 69 80 sig(llm+1)=0. 70 71 DO l = 1, llm 72 dsig(l) = sig(l)-sig(l+1) 73 end DO 74 ELSE 75 if (ok_strato) then 76 if (llm==39) then 77 dsigmin=0.3 78 else if (llm==50) then 79 dsigmin=1. 80 else 81 write(lunout,*) trim(modname), & 82 ' ATTENTION discretisation z a ajuster' 83 dsigmin=1. 84 endif 85 write(lunout,*) trim(modname), & 86 ' Discretisation verticale DSIGMIN=',dsigmin 87 endif 88 81 case("tropo") 89 82 DO l = 1, llm 90 83 x = 2*asin(1.) * (l - 0.5) / (llm + 1) 91 92 IF (ok_strato) THEN 93 dsig(l) =(dsigmin + 7. * SIN(x)**2) & 94 *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 95 ELSE 96 dsig(l) = 1.0 + 7.0 * SIN(x)**2 97 ENDIF 84 dsig(l) = 1.0 + 7.0 * SIN(x)**2 98 85 ENDDO 99 86 dsig = dsig / sum(dsig) … … 102 89 sig(l) = sig(l+1) + dsig(l) 103 90 ENDDO 104 ENDIF 91 case("strato") 92 if (llm==39) then 93 dsigmin=0.3 94 else if (llm==50) then 95 dsigmin=1. 96 else 97 write(lunout,*) trim(modname), ' ATTENTION discretisation z a ajuster' 98 dsigmin=1. 99 endif 100 WRITE(LUNOUT,*) trim(modname), 'Discretisation verticale DSIGMIN=',dsigmin 101 102 DO l = 1, llm 103 x = 2*asin(1.) * (l - 0.5) / (llm + 1) 104 dsig(l) =(dsigmin + 7. * SIN(x)**2) & 105 *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2 106 ENDDO 107 dsig = dsig / sum(dsig) 108 sig(llm+1) = 0. 109 DO l = llm, 1, -1 110 sig(l) = sig(l+1) + dsig(l) 111 ENDDO 112 case("read") 113 call new_unit(unit) 114 open(unit, file="hybrid.txt", status="old", action="read", & 115 position="rewind") 116 read(unit, fmt=*) ! skip title line 117 do l = 1, llm + 1 118 read(unit, fmt=*) sig(l) 119 end do 120 close(unit) 121 case default 122 call abort_gcm("disvert", 'Wrong value for "vert_sampling"', 1) 123 END select 105 124 106 125 DO l=1, llm -
LMDZ5/branches/testing/libf/dyn3d/dynetat0.F
r1421 r1665 119 119 day_ini = tab_cntrl(30) 120 120 itau_dyn = tab_cntrl(31) 121 start_time = tab_cntrl(32) 121 122 c ................................................................. 122 123 c -
LMDZ5/branches/testing/libf/dyn3d/dynredem.F
r1664 r1665 120 120 tab_cntrl(30) = REAL(iday_end) 121 121 tab_cntrl(31) = REAL(itau_dyn + itaufin) 122 c start_time: start_time of simulation (not necessarily 0.) 123 tab_cntrl(32) = start_time 122 124 c 123 125 c ......................................................... -
LMDZ5/branches/testing/libf/dyn3d/etat0_netcdf.F90
r1520 r1665 251 251 !******************************************************************************* 252 252 CALL pression(ip1jmp1, ap, bp, psol, p3d) 253 if ( disvert_type.eq.1) then253 if (pressure_exner) then 254 254 CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y) 255 else ! we assume that we are in the disvert_type==2 case255 else 256 256 CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y) 257 257 endif -
LMDZ5/branches/testing/libf/dyn3d/exner_hyb.F
r1520 r1665 56 56 ! Sanity check 57 57 if (firstcall) then 58 ! check that vertical discretization is compatible59 ! with this routine60 if (disvert_type.ne.1) then61 call abort_gcm(modname,62 & "this routine should only be called if disvert_type==1",42)63 endif64 65 58 ! sanity checks for Shallow Water case (1 vertical layer) 66 59 if (llm.eq.1) then -
LMDZ5/branches/testing/libf/dyn3d/exner_milieu.F
r1520 r1665 53 53 ! Sanity check 54 54 if (firstcall) then 55 ! check that vertical discretization is compatible56 ! with this routine57 if (disvert_type.ne.2) then58 call abort_gcm(modname,59 & "this routine should only be called if disvert_type==2",42)60 endif61 62 55 ! sanity checks for Shallow Water case (1 vertical layer) 63 56 if (llm.eq.1) then -
LMDZ5/branches/testing/libf/dyn3d/gcm.F
r1664 r1665 21 21 ! A nettoyer. On ne veut qu'une ou deux routines d'interface 22 22 ! dynamique -> physique pour l'initialisation 23 ! Ehouarn: for now these only apply to Earth: 24 #ifdef CPP_EARTH 23 #ifdef CPP_PHYS 25 24 USE dimphy 26 25 USE comgeomphy … … 88 87 89 88 REAL zdtvr 90 INTEGER nbetatmoy, nbetatdem,nbetat91 89 92 90 c variables dynamiques … … 181 179 ! A nettoyer. On ne veut qu'une ou deux routines d'interface 182 180 ! dynamique -> physique pour l'initialisation 183 ! Ehouarn : temporarily (?) keep this only for Earth 184 if (planet_type.eq."earth") then 185 #ifdef CPP_EARTH 181 #ifdef CPP_PHYS 186 182 CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) 187 183 call InitComgeomphy 188 184 #endif 189 endif190 185 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 191 186 c----------------------------------------------------------------------- … … 305 300 C on remet le calendrier à zero si demande 306 301 c 302 IF (start_time /= starttime) then 303 WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' 304 &,' fichier restart ne correspond pas à celle lue dans le run.def' 305 IF (raz_date == 1) then 306 WRITE(lunout,*)'Je prends l''heure lue dans run.def' 307 start_time = starttime 308 ELSE 309 WRITE(lunout,*)'Je m''arrete' 310 CALL abort 311 ENDIF 312 ENDIF 307 313 IF (raz_date == 1) THEN 308 314 annee_ref = anneeref … … 373 379 #endif 374 380 375 c nombre d'etats dans les fichiers demarrage et histoire376 nbetatdem = nday / iecri377 nbetatmoy = nday / periodav + 1378 381 379 382 if (iflag_phys.eq.1) then … … 428 431 WRITE(lunout,*) 429 432 . 'GCM: WARNING!!! vitesse verticale nulle dans la physique' 430 ! Earth: 431 if (planet_type.eq."earth") then 432 #ifdef CPP_EARTH 433 ! Physics: 434 #ifdef CPP_PHYS 433 435 CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys , 434 436 , latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp ) 435 437 #endif 436 endif ! of if (planet_type.eq."earth")437 438 call_iniphys=.false. 438 439 ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1)) 439 !#endif440 440 441 441 c numero de stockage pour les fichiers de redemarrage: … … 459 459 #endif 460 460 461 #ifdef CPP_ EARTH461 #ifdef CPP_PHYS 462 462 ! Create start file (startphy.nc) and boundary conditions (limit.nc) 463 463 ! for the Earth verstion -
LMDZ5/branches/testing/libf/dyn3d/guide_mod.F90
r1520 r1665 644 644 ! ----------------------------------------------------------------- 645 645 CALL pression( ip1jmp1, ap, bp, psi, p ) 646 if ( disvert_type==1) then646 if (pressure_exner) then 647 647 CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf) 648 else ! we assume that we are in the disvert_type==2 case648 else 649 649 CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf) 650 650 endif -
LMDZ5/branches/testing/libf/dyn3d/iniacademic.F90
r1664 r1665 222 222 223 223 CALL pression ( ip1jmp1, ap, bp, ps, p ) 224 if ( disvert_type.eq.1) then224 if (pressure_exner) then 225 225 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 226 else if (disvert_type.eq.2) then226 else 227 227 call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf) 228 else229 write(abort_message,*) "Wrong value for disvert_type: ", &230 disvert_type231 call abort_gcm(modname,abort_message,0)232 228 endif 233 229 CALL massdair(p,masse) -
LMDZ5/branches/testing/libf/dyn3d/inidissip.F90
r1502 r1665 28 28 ! Local variables: 29 29 REAL fact,zvert(llm),zz 30 REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm) 30 REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1) 31 real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm) 31 32 REAL ullm,vllm,umin,vmin,zhmin,zhmax 32 REAL zllm ,z1llm33 REAL zllm 33 34 34 35 INTEGER l,ij,idum,ii … … 78 79 DO l = 1,50 79 80 IF(lstardis) THEN 80 CALL divgrad2(1,zh,deltap,niterh, zh)81 CALL divgrad2(1,zh,deltap,niterh,divgra) 81 82 ELSE 82 CALL divgrad (1,zh,niterh, zh)83 CALL divgrad (1,zh,niterh,divgra) 83 84 ENDIF 84 85 85 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 86 87 zllm = ABS( zhmax ) 88 z1llm = 1./zllm 89 DO ij = 1,ip1jmp1 90 zh(ij) = zh(ij)* z1llm 91 ENDDO 86 zllm = ABS(maxval(divgra)) 87 zh = divgra / zllm 92 88 ENDDO 93 89 … … 123 119 !cccc CALL covcont( 1,zu,zv,zu,zv ) 124 120 IF(lstardis) THEN 125 CALL gradiv2( 1,zu,zv,nitergdiv, zu,zv)121 CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy ) 126 122 ELSE 127 CALL gradiv ( 1,zu,zv,nitergdiv, zu,zv)123 CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy ) 128 124 ENDIF 129 125 ELSE 130 126 IF(lstardis) THEN 131 CALL nxgraro2( 1,zu,zv,nitergrot, zu,zv)127 CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy ) 132 128 ELSE 133 CALL nxgrarot( 1,zu,zv,nitergrot, zu,zv)129 CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy ) 134 130 ENDIF 135 131 ENDIF 136 132 137 CALL minmax(iip1*jjp1,zu,umin,ullm ) 138 CALL minmax(iip1*jjm, zv,vmin,vllm ) 139 140 ullm = ABS ( ullm ) 141 vllm = ABS ( vllm ) 142 143 zllm = MAX( ullm,vllm ) 144 z1llm = 1./ zllm 145 DO ij = 1, ip1jmp1 146 zu(ij) = zu(ij)* z1llm 147 ENDDO 148 DO ij = 1, ip1jm 149 zv(ij) = zv(ij)* z1llm 150 ENDDO 133 zllm = max(abs(maxval(gx)), abs(maxval(gy))) 134 zu = gx / zllm 135 zv = gy / zllm 151 136 end DO 152 137 -
LMDZ5/branches/testing/libf/dyn3d/inigrads.F
r524 r1665 9 9 implicit none 10 10 11 integer if,im,jm,lm,i,j,l ,lnblnk11 integer if,im,jm,lm,i,j,l 12 12 real x(im),y(jm),z(lm),fx,fy,fz,dt 13 13 real xmin,xmax,ymin,ymax … … 40 40 ivar(if)=0 41 41 42 fichier(if)= file(1:lnblnk(file))42 fichier(if)=trim(file) 43 43 44 44 firsttime(if)=.true. … … 70 70 71 71 print*,4*(ifd(if)-iid(if))*(jfd(if)-jid(if)) 72 print*, file(1:lnblnk(file))//'.dat'72 print*,trim(file)//'.dat' 73 73 74 OPEN (unit(if)+1,FILE= file(1:lnblnk(file))//'.dat'74 OPEN (unit(if)+1,FILE=trim(file)//'.dat' 75 75 s ,FORM='unformatted', 76 76 s ACCESS='direct' -
LMDZ5/branches/testing/libf/dyn3d/integrd.F
r1550 r1665 4 4 SUBROUTINE integrd 5 5 $ ( nq,vcovm1,ucovm1,tetam1,psm1,massem1, 6 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis,finvmaold ) 6 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis !,finvmaold 7 & ) 7 8 8 9 use control_mod, only : planet_type … … 34 35 #include "temps.h" 35 36 #include "serre.h" 37 #include "iniprint.h" 36 38 37 39 c Arguments: 38 40 c ---------- 39 41 40 INTEGER nq 41 42 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm) 43 REAL q(ip1jmp1,llm,nq) 44 REAL ps(ip1jmp1),masse(ip1jmp1,llm),phis(ip1jmp1) 45 46 REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm) 47 REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1),massem1(ip1jmp1,llm) 48 49 REAL dv(ip1jm,llm),du(ip1jmp1,llm) 50 REAL dteta(ip1jmp1,llm),dp(ip1jmp1) 51 REAL dq(ip1jmp1,llm,nq), finvmaold(ip1jmp1,llm) 42 integer,intent(in) :: nq ! number of tracers to handle in this routine 43 real,intent(inout) :: vcov(ip1jm,llm) ! covariant meridional wind 44 real,intent(inout) :: ucov(ip1jmp1,llm) ! covariant zonal wind 45 real,intent(inout) :: teta(ip1jmp1,llm) ! potential temperature 46 real,intent(inout) :: q(ip1jmp1,llm,nq) ! advected tracers 47 real,intent(inout) :: ps(ip1jmp1) ! surface pressure 48 real,intent(inout) :: masse(ip1jmp1,llm) ! atmospheric mass 49 real,intent(in) :: phis(ip1jmp1) ! ground geopotential !!! unused 50 ! values at previous time step 51 real,intent(inout) :: vcovm1(ip1jm,llm) 52 real,intent(inout) :: ucovm1(ip1jmp1,llm) 53 real,intent(inout) :: tetam1(ip1jmp1,llm) 54 real,intent(inout) :: psm1(ip1jmp1) 55 real,intent(inout) :: massem1(ip1jmp1,llm) 56 ! the tendencies to add 57 real,intent(in) :: dv(ip1jm,llm) 58 real,intent(in) :: du(ip1jmp1,llm) 59 real,intent(in) :: dteta(ip1jmp1,llm) 60 real,intent(in) :: dp(ip1jmp1) 61 real,intent(in) :: dq(ip1jmp1,llm,nq) !!! unused 62 ! real,intent(out) :: finvmaold(ip1jmp1,llm) !!! unused 52 63 53 64 c Local: … … 55 66 56 67 REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1) 57 REAL massescr( ip1jmp1,llm ), finvmasse(ip1jmp1,llm) 68 REAL massescr( ip1jmp1,llm ) 69 ! REAL finvmasse(ip1jmp1,llm) 58 70 REAL p(ip1jmp1,llmp1) 59 71 REAL tpn,tps,tppn(iim),tpps(iim) … … 61 73 REAL deltap( ip1jmp1,llm ) 62 74 63 INTEGER l,ij,iq 75 INTEGER l,ij,iq,i,j 64 76 65 77 REAL SSUM … … 88 100 DO ij = 1,ip1jmp1 89 101 IF( ps(ij).LT.0. ) THEN 90 PRINT*,' Au point ij = ',ij, ' , pression sol neg. ', ps(ij) 91 print *, ' dans integrd' 92 stop 1 102 write(lunout,*) "integrd: negative surface pressure ",ps(ij) 103 write(lunout,*) " at node ij =", ij 104 ! since ij=j+(i-1)*jjp1 , we have 105 j=modulo(ij,jjp1) 106 i=1+(ij-j)/jjp1 107 write(lunout,*) " lon = ",rlonv(i)*180./pi, " deg", 108 & " lat = ",rlatu(j)*180./pi, " deg" 109 stop 93 110 ENDIF 94 111 ENDDO … … 110 127 CALL massdair ( p , masse ) 111 128 112 CALL SCOPY( ijp1llm , masse, 1, finvmasse, 1 ) 113 CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1 ) 129 ! Ehouarn : we don't use/need finvmaold and finvmasse, 130 ! so might as well not compute them 131 ! CALL SCOPY( ijp1llm , masse, 1, finvmasse, 1 ) 132 ! CALL filtreg( finvmasse, jjp1, llm, -2, 2, .TRUE., 1 ) 114 133 c 115 134 … … 218 237 ENDDO 219 238 220 221 CALL SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )239 ! Ehouarn: forget about finvmaold 240 ! CALL SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 ) 222 241 223 242 endif ! of if (planet_type.eq."earth") -
LMDZ5/branches/testing/libf/dyn3d/leapfrog.F
r1529 r1665 72 72 73 73 real zqmin,zqmax 74 INTEGER nbetatmoy, nbetatdem,nbetat75 74 76 75 c variables dynamiques … … 118 117 119 118 REAL SSUM 120 REAL time_0 , finvmaold(ip1jmp1,llm) 119 REAL time_0 120 ! REAL finvmaold(ip1jmp1,llm) 121 121 122 122 cym LOGICAL lafin … … 212 212 dq(:,:,:)=0. 213 213 CALL pression ( ip1jmp1, ap, bp, ps, p ) 214 if ( disvert_type==1) then214 if (pressure_exner) then 215 215 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 216 else ! we assume that we are in the disvert_type==2 case216 else 217 217 CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) 218 218 endif … … 222 222 c ---------------------------------- 223 223 224 1 CONTINUE 225 226 jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec) 227 jH_cur = jH_ref + & 228 & (itau * dtvr / daysec - int(itau * dtvr / daysec)) 224 1 CONTINUE ! Matsuno Forward step begins here 225 226 jD_cur = jD_ref + day_ini - day_ref + & 227 & itau/day_step 228 jH_cur = jH_ref + start_time + & 229 & mod(itau,day_step)/float(day_step) 230 jD_cur = jD_cur + int(jH_cur) 231 jH_cur = jH_cur - int(jH_cur) 229 232 230 233 … … 255 258 256 259 c ... P.Le Van .26/04/94 .... 257 258 CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )259 CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )260 261 2 CONTINUE 260 ! Ehouarn: finvmaold is actually not used 261 ! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 ) 262 ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 ) 263 264 2 CONTINUE ! Matsuno backward or leapfrog step begins here 262 265 263 266 c----------------------------------------------------------------------- … … 302 305 c -------------------------------- 303 306 307 ! compute geopotential phi() 304 308 CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi ) 305 309 … … 315 319 316 320 IF( forward. OR . leapf ) THEN 317 321 ! Ehouarn: NB: at this point p with ps are not synchronized 322 ! (whereas mass and ps are...) 318 323 CALL caladvtrac(q,pbaru,pbarv, 319 324 * p, masse, dq, teta, … … 340 345 341 346 CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 , 342 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,343 $ finvmaold )347 $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ) 348 ! $ finvmaold ) 344 349 345 350 … … 364 369 365 370 CALL pression ( ip1jmp1, ap, bp, ps, p ) 366 if ( disvert_type==1) then371 if (pressure_exner) then 367 372 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf ) 368 else ! we assume that we are in the disvert_type==2 case373 else 369 374 CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) 370 375 endif … … 372 377 ! rdaym_ini = itau * dtvr / daysec 373 378 ! rdayvrai = rdaym_ini + day_ini 374 jD_cur = jD_ref + day_ini - day_ref 375 $ + int (itau * dtvr / daysec) 376 jH_cur = jH_ref + & 377 & (itau * dtvr / daysec - int(itau * dtvr / daysec)) 379 ! jD_cur = jD_ref + day_ini - day_ref 380 ! $ + int (itau * dtvr / daysec) 381 ! jH_cur = jH_ref + & 382 ! & (itau * dtvr / daysec - int(itau * dtvr / daysec)) 383 jD_cur = jD_ref + day_ini - day_ref + & 384 & itau/day_step 385 jH_cur = jH_ref + start_time + & 386 & mod(itau,day_step)/float(day_step) 387 jD_cur = jD_cur + int(jH_cur) 388 jH_cur = jH_cur - int(jH_cur) 378 389 ! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur 379 390 ! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes) … … 394 405 ! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)! 395 406 IF (planet_type.eq."earth") THEN 407 #ifdef CPP_EARTH 396 408 CALL diagedyn(ztit,2,1,1,dtphys 397 409 & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2)) 410 #endif 398 411 ENDIF 399 412 ENDIF ! of IF (ip_ebil_dyn.ge.1 ) … … 472 485 473 486 CALL pression ( ip1jmp1, ap, bp, ps, p ) 474 if ( disvert_type==1) then487 if (pressure_exner) then 475 488 CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) 476 else ! we assume that we are in the disvert_type==2 case489 else 477 490 CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf ) 478 491 endif … … 529 542 ENDDO 530 543 531 DO ij = 1,iim 532 tppn(ij) = aire( ij ) * ps ( ij ) 533 tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm) 534 ENDDO 535 tpn = SSUM(iim,tppn,1)/apoln 536 tps = SSUM(iim,tpps,1)/apols 537 538 DO ij = 1, iip1 539 ps( ij ) = tpn 540 ps(ij+ip1jm) = tps 541 ENDDO 542 544 if (1 == 0) then 545 !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics 546 !!! 2) should probably not be here anyway 547 !!! but are kept for those who would want to revert to previous behaviour 548 DO ij = 1,iim 549 tppn(ij) = aire( ij ) * ps ( ij ) 550 tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm) 551 ENDDO 552 tpn = SSUM(iim,tppn,1)/apoln 553 tps = SSUM(iim,tpps,1)/apols 554 555 DO ij = 1, iip1 556 ps( ij ) = tpn 557 ps(ij+ip1jm) = tps 558 ENDDO 559 endif ! of if (1 == 0) 543 560 544 561 END IF ! of IF(apdiss) … … 622 639 ! Ehouarn: output only during LF or Backward Matsuno 623 640 if (leapf.or.(.not.leapf.and.(.not.forward))) then 624 nbetat = nbetatdem625 641 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) 626 642 unat=0. … … 652 668 ! if (planet_type.eq."earth") then 653 669 ! Write an Earth-format restart file 654 CALL dynredem1("restart.nc", 0.0,670 CALL dynredem1("restart.nc",start_time, 655 671 & vcov,ucov,teta,q,masse,ps) 656 672 ! endif ! of if (planet_type.eq."earth") 657 673 658 674 CLOSE(99) 675 !!! Ehouarn: Why not stop here and now? 659 676 ENDIF ! of IF (itau.EQ.itaufin) 660 677 … … 740 757 IF(MOD(itau,iecri ).EQ.0) THEN 741 758 c IF(MOD(itau,iecri*day_step).EQ.0) THEN 742 nbetat = nbetatdem743 759 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) 744 760 unat=0. … … 763 779 IF(itau.EQ.itaufin) THEN 764 780 ! if (planet_type.eq."earth") then 765 CALL dynredem1("restart.nc", 0.0,781 CALL dynredem1("restart.nc",start_time, 766 782 & vcov,ucov,teta,q,masse,ps) 767 783 ! endif ! of if (planet_type.eq."earth") -
LMDZ5/branches/testing/libf/dyn3d/temps.h
r1279 r1665 14 14 15 15 COMMON/temps/itaufin, dt, day_ini, day_end, annee_ref, day_ref, & 16 & itau_dyn, itau_phy, jD_ref, jH_ref, calend 16 & itau_dyn, itau_phy, jD_ref, jH_ref, calend, & 17 & start_time 18 17 19 18 20 INTEGER itaufin 19 21 INTEGER itau_dyn, itau_phy 20 22 INTEGER day_ini, day_end, annee_ref, day_ref 21 REAL dt, jD_ref, jH_ref 23 REAL dt, jD_ref, jH_ref, start_time 22 24 CHARACTER (len=10) :: calend 23 25 -
LMDZ5/branches/testing/libf/dyn3d/wrgrads.F
r1025 r1665 26 26 c local 27 27 28 integer im,jm,lm,i,j,l, lnblnk,iv,iii,iji,iif,ijf28 integer im,jm,lm,i,j,l,iv,iii,iji,iif,ijf 29 29 30 30 logical writectl … … 59 59 nvar(if)=ivar(if) 60 60 var(ivar(if),if)=name 61 tvar(ivar(if),if)=t itlevar(1:lnblnk(titlevar))61 tvar(ivar(if),if)=trim(titlevar) 62 62 nld(ivar(if),if)=nl 63 63 c print*,'initialisation ecriture de ',var(ivar(if),if) … … 101 101 file=fichier(if) 102 102 c WARNING! on reecrase le fichier .ctl a chaque ecriture 103 open(unit(if),file= file(1:lnblnk(file))//'.ctl'103 open(unit(if),file=trim(file)//'.ctl' 104 104 & ,form='formatted',status='unknown') 105 105 write(unit(if),'(a5,1x,a40)') 106 & 'DSET ','^'// file(1:lnblnk(file))//'.dat'106 & 'DSET ','^'//trim(file)//'.dat' 107 107 108 108 write(unit(if),'(a12)') 'UNDEF 1.0E30'
Note: See TracChangeset
for help on using the changeset viewer.
