Changeset 232


Ignore:
Timestamp:
Jun 20, 2001, 3:53:15 PM (23 years ago)
Author:
lmdzadmin
Message:

Merge par rapport a la branche principale
LF

Location:
LMDZ.3.3/branches/rel-LF/libf/dyn3d
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/calfis.F

    r177 r232  
    319319c
    320320      DO l=1,llm
    321         pvervel(1,l)=pw(1,1,l)/apoln
     321        pvervel(1,l)=pw(1,1,l) * g /apoln
    322322        ig0=2
    323323       DO j=2,jjm
    324324           DO i = 1, iim
    325               pvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
     325              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
    326326              ig0 = ig0 + 1
    327327           ENDDO
    328328       ENDDO
    329         pvervel(ig0,l)=pw(1,jjp1,l)/apols
     329        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
    330330      ENDDO
    331331

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/comvert.h

    r2 r232  
    33
    44      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm) ,
    5      ,               pa,preff,nivsigs(llm),nivsig(llmp1)
     5     ,               pa,preff,nivsigs(llm),nivsig(llm+1)
    66
    77      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F

    r177 r232  
    305305        ENDIF
    306306      END DO
    307 
    308307c
    309308c
    … …  
    338337        STOP
    339338      ENDIF
     339#ifdef NC_DOUBLE
     340      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
     341#else
    340342      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
     343#endif
    341344      if (ierr.ne.0) then
    342345        print *, NF_STRERROR(ierr)
    … …  
    354357        STOP
    355358      ENDIF
     359#ifdef NC_DOUBLE
     360      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
     361#else
    356362      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
     363#endif
    357364      if (ierr.ne.0) then
    358365        print *, NF_STRERROR(ierr)
    … …  
    370377        STOP
    371378      ENDIF
     379#ifdef NC_DOUBLE
     380      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
     381#else
    372382      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
     383#endif
    373384      if (ierr.ne.0) then
    374385        print *, NF_STRERROR(ierr)
    … …  
    387398         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
    388399         print*,dimfirst,dimlast
     400#ifdef NC_DOUBLE
     401         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
     402#else
    389403         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
     404#endif
    390405         if (ierr.ne.0) then
    391406           print *, NF_STRERROR(ierr)
    … …  
    471486        STOP
    472487      ENDIF
     488#ifdef NC_DOUBLE
     489      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
     490#else
    473491      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
     492#endif
    474493      if (ierr.ne.0) then
    475494        print *, NF_STRERROR(ierr)
    … …  
    487506        STOP
    488507      ENDIF
     508#ifdef NC_DOUBLE
     509      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
     510#else
    489511      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
     512#endif
    490513      if (ierr.ne.0) then
    491514        print *, NF_STRERROR(ierr)
    … …  
    503526        STOP
    504527      ENDIF
     528#ifdef NC_DOUBLE
     529      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
     530#else
    505531      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
     532#endif
    506533      if (ierr.ne.0) then
    507534        print *, NF_STRERROR(ierr)
    … …  
    519546c
    520547         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
     548#ifdef NC_DOUBLE
     549         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
     550#else
    521551         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
     552#endif
    522553         if (ierr.ne.0) then
    523554           print *, NF_STRERROR(ierr)
    … …  
    698729        STOP
    699730      ENDIF
     731#ifdef NC_DOUBLE
     732      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
     733#else
    700734      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
     735#endif
     736
    701737      if (ierr.ne.0) then
    702738        print *, NF_STRERROR(ierr)
    … …  
    714750        STOP
    715751      ENDIF
     752#ifdef NC_DOUBLE
     753      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
     754#else
    716755      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
     756#endif
    717757      if (ierr.ne.0) then
    718758        print *, NF_STRERROR(ierr)
    … …  
    730770        STOP
    731771      ENDIF
     772#ifdef NC_DOUBLE
     773      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
     774#else
    732775      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
     776#endif
    733777      if (ierr.ne.0) then
    734778        print *, NF_STRERROR(ierr)
    … …  
    746790c
    747791         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
     792#ifdef NC_DOUBLE
     793         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
     794#else
    748795         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
     796#endif
    749797         if (ierr.ne.0) then
    750798           print *, NF_STRERROR(ierr)
    … …  
    825873        STOP
    826874      ENDIF
     875#ifdef NC_DOUBLE
     876      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
     877#else
    827878      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
     879#endif
    828880      if (ierr.ne.0) then
    829881        print *, NF_STRERROR(ierr)
    … …  
    841893        STOP
    842894      ENDIF
     895#ifdef NC_DOUBLE
     896      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
     897#else
    843898      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
     899#endif
    844900      if (ierr.ne.0) then
    845901        print *, NF_STRERROR(ierr)
    … …  
    857913        STOP
    858914      ENDIF
     915#ifdef NC_DOUBLE
     916      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
     917#else
    859918      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
     919#endif
    860920      if (ierr.ne.0) then
    861921        print *, NF_STRERROR(ierr)
    … …  
    873933c
    874934         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
     935#ifdef NC_DOUBLE
     936         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
     937#else
    875938         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
     939#endif
    876940         if (ierr.ne.0) then
    877941           print *, NF_STRERROR(ierr)

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/defrun_new.F

    r177 r232  
     1c
     2c $Header$
     3c
    14      SUBROUTINE defrun_new( tapedef, etatinit, clesphy0 )
    25c
    … …  
    3639      INTEGER   tapeout
    3740      REAL clonn,clatt,grossismxx,grossismyy
    38       REAL dzoomxx,dzoomyy
     41      REAL dzoomxx,dzoomyy,tauxx,tauyy
    3942      LOGICAL  fxyhypbb, ysinuss
    4043      INTEGER i
    … …  
    258261      WRITE(tapeout,*)    clonn
    259262      IF( ABS(clon - clonn).GE. 0.001 )  THEN
    260         PRINT *,' La valeur de clon passee par run.def est differente de
    261      *  celle lue sur le fichier  start '
     263       WRITE(tapeout,*) ' La valeur de clon passee par run.def est diffe
     264     *rente de  celle lue sur le fichier  start '
    262265        STOP
    263266      ENDIF
    264 c
    265267c
    266268      READ (tapedef,9001) ch1,ch4
    … …  
    270272
    271273      IF( ABS(clat - clatt).GE. 0.001 )  THEN
    272         PRINT *,' La valeur de clat passee par run.def est differente de
    273      *  celle lue sur le fichier  start '
     274       WRITE(tapeout,*) ' La valeur de clat passee par run.def est diffe
     275     *rente de  celle lue sur le fichier  start '
    274276        STOP
    275277      ENDIF
    … …  
    281283
    282284      IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
    283         PRINT *,' La valeur de grossismx passee par run.def est differente 
    284      * de celle lue sur le fichier  start '
     285       WRITE(tapeout,*) ' La valeur de grossismx passee par run.def est
     286     , differente de celle lue sur le fichier  start '
    285287        STOP
    286288      ENDIF
    … …  
    292294
    293295      IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
    294         PRINT *,' La valeur de grossismy passee par run.def est differen
    295      * te de celle lue sur le fichier  start '
     296       WRITE(tapeout,*) ' La valeur de grossismy passee par run.def est
     297     , differente de celle lue sur le fichier  start '
    296298        STOP
    297299      ENDIF
    298300     
    299301      IF( grossismx.LT.1. )  THEN
    300         PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
     302        WRITE(tapeout,*) ' ***  ATTENTION !! grossismx < 1 .   *** '
    301303         STOP
    302304      ELSE
    … …  
    306308
    307309      IF( grossismy.LT.1. )  THEN
    308         PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
     310        WRITE(tapeout,*) ' ***  ATTENTION !! grossismy < 1 .   *** '
    309311         STOP
    310312      ELSE
    … …  
    312314      ENDIF
    313315
    314       PRINT *,' alphax alphay defrun ',alphax,alphay
    315316c
    316317c    alphax et alphay sont les anciennes formulat. des grossissements
    … …  
    324325      IF( .NOT.fxyhypb )  THEN
    325326           IF( fxyhypbb )     THEN
    326               PRINT *,' ********  PBS DANS  DEFRUN  ******** '
    327               PRINT *,' *** fxyhypb lu sur le fichier start est F ',
    328      *       'alors  qu il est  T  sur  run.def  ***'
     327            WRITE(tapeout,*) ' ********  PBS DANS  DEFRUN  ******** '
     328            WRITE(tapeout,*)' *** fxyhypb lu sur le fichier start est F'
     329     *,      '                   alors  qu il est  T  sur  run.def  ***'
    329330              STOP
    330331           ENDIF
    331332      ELSE
    332333           IF( .NOT.fxyhypbb )   THEN
    333               PRINT *,' ********  PBS DANS  DEFRUN  ******** '
    334               PRINT *,' ***  fxyhypb lu sur le fichier start est T ',
    335      *        'alors  qu il est  F  sur  run.def  ****  '
     334            WRITE(tapeout,*) ' ********  PBS DANS  DEFRUN  ******** '
     335            WRITE(tapeout,*)' *** fxyhypb lu sur le fichier start est t'
     336     *,      '                   alors  qu il est  F  sur  run.def  ***'
    336337              STOP
    337338           ENDIF
    … …  
    343344      WRITE(tapeout,*)    dzoomxx
    344345
    345       IF( fxyhypb )  THEN
    346        IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
    347         PRINT *,' La valeur de dzoomx passee par run.def est differente
    348      *  de celle lue sur le fichier  start '
    349         STOP
    350        ENDIF
    351       ENDIF
    352 
    353346      READ (tapedef,9001) ch1,ch4
    354347      READ (tapedef,*)    dzoomyy
    … …  
    356349      WRITE(tapeout,*)    dzoomyy
    357350
     351      READ (tapedef,9001) ch1,ch4
     352      READ (tapedef,*)    tauxx
     353      WRITE(tapeout,9001) ch1,'taux'
     354      WRITE(tapeout,*)    tauxx
     355
     356      READ (tapedef,9001) ch1,ch4
     357      READ (tapedef,*)    tauyy
     358      WRITE(tapeout,9001) ch1,'tauy'
     359      WRITE(tapeout,*)    tauyy
     360
    358361      IF( fxyhypb )  THEN
     362
     363       IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
     364        WRITE(tapeout,*)' La valeur de dzoomx passee par run.def est dif
     365     *ferente de celle lue sur le fichier  start '
     366        CALL ABORT
     367       ENDIF
     368
    359369       IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
    360         PRINT *,' La valeur de dzoomy passee par run.def est differente
    361      * de celle lue sur le fichier  start '
    362         STOP
     370        WRITE(tapeout,*)' La valeur de dzoomy passee par run.def est dif
     371     *ferente de celle lue sur le fichier  start '
     372        CALL ABORT
    363373       ENDIF
     374
     375       IF( ABS(taux - tauxx).GE. 0.001 )  THEN
     376        WRITE(6,*)' La valeur de taux passee par run.def est differente
     377     *  de celle lue sur le fichier  start '
     378        CALL ABORT
     379       ENDIF
     380
     381       IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
     382        WRITE(6,*)' La valeur de tauy passee par run.def est differente
     383     *  de celle lue sur le fichier  start '
     384        CALL ABORT
     385       ENDIF
     386
    364387      ENDIF
    365388     
    … …  
    374397        IF( .NOT.ysinus )  THEN
    375398           IF( ysinuss )     THEN
    376               PRINT *,' ********  PBS DANS  DEFRUN  ******** '
    377               PRINT *,' *** ysinus lu sur le fichier start est F ',
    378      *       'alors  qu il est  T  sur  run.def  ***'
     399              WRITE(6,*) ' ********  PBS DANS  DEFRUN  ******** '
     400              WRITE(tapeout,*)'** ysinus lu sur le fichier start est F',
     401     *       ' alors  qu il est  T  sur  run.def  ***'
    379402              STOP
    380403           ENDIF
    381404        ELSE
    382405           IF( .NOT.ysinuss )   THEN
    383               PRINT *,' ********  PBS DANS  DEFRUN  ******** '
    384               PRINT *,' ***  ysinus lu sur le fichier start est T ',
    385      *        'alors  qu il est  F  sur  run.def  ****  '
     406              WRITE(6,*) ' ********  PBS DANS  DEFRUN  ******** '
     407              WRITE(tapeout,*)'** ysinus lu sur le fichier start est T',
     408     *       ' alors  qu il est  F  sur  run.def  ***'
    386409              STOP
    387410           ENDIF
    … …  
    389412      ENDIF
    390413c
    391       close(tapedef)
     414      WRITE(6,*) ' alphax alphay defrun ',alphax,alphay
     415
     416      CLOSE(tapedef)
     417
    392418      RETURN
    393419c   ...............................................
    … …  
    416442
    417443      IF( grossismx.LT.1. )  THEN
    418         PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
     444        WRITE(tapeout,*) '***  ATTENTION !! grossismx < 1 .   *** '
    419445         STOP
    420446      ELSE
    … …  
    423449
    424450      IF( grossismy.LT.1. )  THEN
    425         PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
     451        WRITE(tapeout,*) ' ***  ATTENTION !! grossismy < 1 .   *** '
    426452         STOP
    427453      ELSE
    … …  
    429455      ENDIF
    430456
    431       PRINT *,' alphax alphay defrun ',alphax,alphay
    432 
    433457c
    434458      READ (tapedef,9001) ch1,ch4
    … …  
    446470      WRITE(tapeout,9001) ch1,'dzoomy'
    447471      WRITE(tapeout,*)    dzoomy
    448 c
     472
     473      READ (tapedef,9001) ch1,ch4
     474      READ (tapedef,*)    taux
     475      WRITE(tapeout,9001) ch1,'taux'
     476      WRITE(tapeout,*)    taux
     477c
     478      READ (tapedef,9001) ch1,ch4
     479      READ (tapedef,*)    tauy
     480      WRITE(tapeout,9001) ch1,'tauy'
     481      WRITE(tapeout,*)    tauy
     482
    449483      READ (tapedef,9001) ch1,ch4
    450484      READ (tapedef,*)    ysinus
    … …  
    452486      WRITE(tapeout,*)    ysinus
    453487       
     488      WRITE(tapeout,*) ' alphax alphay defrun ',alphax,alphay
    454489c
    4554909000  FORMAT(3(/,a72))
    4564919001  FORMAT(/,a72,/,a12)
    457492cc
    458       close(tapedef)
     493      CLOSE(tapedef)
     494
    459495      RETURN
    460496      END

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/dynetat0.F

    r2 r232  
     1c
     2c $Header$
     3c
    14      SUBROUTINE dynetat0(fichnom,nq,vcov,ucov,
    25     .                    teta,q,masse,ps,phis,time)
    … …  
    105108        dzoomx   = tab_cntrl(25)
    106109        dzoomy   = tab_cntrl(26)
     110        taux     = tab_cntrl(28)
     111        tauy     = tab_cntrl(29)
    107112      ELSE
    108113        fxyhypb = . FALSE .

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/dynredem.F

    r179 r232  
    9494       tab_cntrl(25) = dzoomx
    9595       tab_cntrl(26) = dzoomy
     96       tab_cntrl(27) = 0.
     97       tab_cntrl(28) = taux
     98       tab_cntrl(29) = tauy
    9699      ELSE
    97100       tab_cntrl(24) = 0.
    … …  
    99102       tab_cntrl(26) = dzoomy
    100103       tab_cntrl(27) = 0.
    101        IF( ysinus )  tab_cntrl(27) = 1.
     104       tab_cntrl(28) = 0.
     105       tab_cntrl(29) = 0.
     106       IF( ysinus )  tab_cntrl(27) = 1.
    102107      ENDIF
    103108c

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F

    r227 r232  
    274274      q3d(:,:,:,:) = 0.0
    275275      qd(:,:,:) = 0.0
     276      q3d(:,:,:,:) = 0.0
    276277      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
    277278     .                           maxval(qsat(:,:,:))

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/fluxstokenc.F

    r79 r232  
    11      SUBROUTINE fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
    2      . time_step,itau, fluxid, fluxvid,fluxdid )
     2     . time_step,itau )
    33
    44       USE IOIPSL
    … …  
    3030
    3131      REAL pbarvst(iip1,jjp1,llm),zistdyn
     32        real dtcum
    3233
    33 
    34       INTEGER iadvtr
    35       integer ndex1d(1)
    36       integer ndex2d(ip1jmp1)
    37       integer ndex3dv(ip1jm*llm),ndex3d(ip1jmp1*llm)
     34      INTEGER iadvtr,ndex(1)
    3835      integer nscal
    3936      real tst(1),ist(1),istp(1)
    … …  
    5350     .  time_step,istdyn* time_step,istdyn* time_step,
    5451     . nqmx, fluxid,fluxvid,fluxdid)
     52       
     53        ndex(1) = 0
     54        call histwrite(fluxid, 'phis', 1, phis, iip1*jjp1, ndex)
     55        call histwrite(fluxid, 'aire', 1, aire, iip1*jjp1, ndex)
     56       
     57        ndex(1) = 0
     58        nscal = 1
     59        tst(1) = time_step
     60        call histwrite(fluxdid, 'dtvr', 1, tst, nscal, ndex)
     61        ist(1)=istdyn
     62        call histwrite(fluxdid, 'istdyn', 1, ist, nscal, ndex)
     63        istp(1)= istphy
     64        call histwrite(fluxdid, 'istphy', 1, istp, nscal, ndex)
     65       
    5566        first = .false.
    5667
    5768      endif
    58 
    59 
    60       ndex1d = 0
    61       ndex2d = 0
    62       ndex3dv = 0
    63       ndex3d = 0
    6469
    6570
    … …  
    9398c   Test pour savoir si on advecte a ce pas de temps
    9499      IF ( iadvtr.EQ.istdyn ) THEN
    95 
    96100c    normalisation
    97101      DO l=1,llm
    … …  
    124128
    125129         iadvtr=0
    126 
    127 c     write(*,*)'histwrite phis'
    128       call histwrite(fluxid, 'phis', 1, phis, iip1*jjp1, ndex2d)
    129 c     write(*,*)'histwrite aire'
    130       call histwrite(fluxid, 'aire', 1, aire, iip1*jjp1, ndex2d)
     130        Print*,'ITAU auqel on stoke les fluxmasses',itau
    131131       
    132       nscal = 1
    133       tst(1) = time_step
    134 c     write(*,*)'histwrite dtvr'
    135       call histwrite(fluxdid, 'dtvr', 1, tst, nscal, ndex1d)
    136       ist(1)=istdyn
    137 c     write(*,*)'histwrite istdyn'
    138       call histwrite(fluxdid, 'istdyn', 1, ist, nscal, ndex1d)
    139       istp(1)= istphy
    140 c     write(*,*)'histwrite istphy'
    141       call histwrite(fluxdid, 'istphy', 1, istp, nscal, ndex1d)
    142 
    143 c       write(*,*)'histwrite masse'     
    144         call histwrite(fluxid, 'masse', itau, masse,
    145      .               iip1*jjp1*llm, ndex3d)
     132        call histwrite(fluxid, 'masse', itau, massem,
     133     .               iip1*jjp1*llm, ndex)
    146134       
    147 c       write(*,*)'histwrite pbaru'     
    148135        call histwrite(fluxid, 'pbaru', itau, pbarug,
    149      .               iip1*jjp1*llm, ndex3d)
     136     .               iip1*jjp1*llm, ndex)
    150137       
    151 c       write(*,*)'histwrite pbarv'     
    152         call histwrite(fluxvid, 'pbarv', itau, pbarvst,
    153      .               iim*jjp1*llm, ndex3dv)
     138        call histwrite(fluxvid, 'pbarv', itau, pbarvg,
     139     .               iip1*jjm*llm, ndex)
    154140       
    155 c       write(*,*)'histwrite w'
    156141        call histwrite(fluxid, 'w' ,itau, wg,
    157      .             iip1*jjp1*llm, ndex3d)
     142     .             iip1*jjp1*llm, ndex)
    158143       
    159 c       write(*,*)'histwrite teta'     
    160144        call histwrite(fluxid, 'teta' ,itau, tetac,
    161      .             iip1*jjp1*llm, ndex3d)
     145     .             iip1*jjp1*llm, ndex)
    162146       
    163 c       write(*,*)'histwrite phi'       
    164147        call histwrite(fluxid, 'phi' ,itau, phic,
    165      .             iip1*jjp1*llm, ndex3d)
     148     .             iip1*jjp1*llm, ndex)
    166149       
    167150C

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxhyp.F

    r2 r232  
    11       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
    2      , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025)
     2     , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,
     3     , champmin,champmax                                               )
    34
    45c      Auteur :  P. Le Van
    … …  
    1112c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
    1213c     dzoom  etant  la distance totale de la zone du zoom
    13 c     tau  la transition ,   normalement  = 1  .
    14 c
     14c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom
     15c
     16c    On doit avoir grossism x dzoom <  pi ( radians )   , en longitude.
     17c   ********************************************************************
     18
    1519
    1620       INTEGER nmax, nmax2
    17        PARAMETER (  nmax = 50000, nmax2 = 2*nmax )
     21       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
    1822
    1923#include "dimensions.h"
    … …  
    2327c
    2428       REAL xzoomdeg,dzoom,tau,grossism
     29
     30c    ......   arguments  de  sortie  ......
     31
    2532       REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
    2633     ,  rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
    2734
    28 c    ......   arguments  de  sortie  ......
    29 c
    30        REAL xlon(iip1),xprimm(iip1),xuv
    31        REAL xtild(0:nmax2)
    32        REAL fhyp(0:nmax),ffdx(0:nmax),beta,Xprimt(0:nmax2)
    33        REAL Xf(0:nmax2),xxpr(0:nmax2)
    34        REAL xvrai(iip1),xxprim(iip1)
    35        REAL pi,depi,epsilon,xzoom
     35c     .... variables locales  ....
     36c
     37       REAL*8 xlon(iip1),xprimm(iip1),xuv
     38       REAL*8 xtild(0:nmax2)
     39       REAL*8 fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)
     40       REAL*8 Xf(0:nmax2),xxpr(0:nmax2)
     41       REAL*8 xvrai(iip1),xxprim(iip1)
     42       REAL*8 pi,depi,epsilon,xzoom,fa,fb
     43       REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2
    3644       INTEGER i,it,ik,iter,ii,idif
    37        REAL xi,xo1,xint,xmoy,xlon2,fxm,Xprimin
    38        REAL champmin,champmax
     45       REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
     46       REAL*8 champmin,champmax
    3947       INTEGER is2
    4048       SAVE is2
    41        REAL dlon1(iip1),dlon2(iip1),dlon3(iip1)
     49
     50       REAL*8 heavyside
     51       EXTERNAL coefpoly,heavyside
    4252
    4353       pi       = 2. * ASIN(1.)
    4454       depi     = 2. * pi
    45        epsilon  = 1.e-6
     55       epsilon  = 1.e-3
    4656       xzoom    = xzoomdeg * pi/180.
    4757
    48 
     58       IF( dzoom.LT.1.)  THEN
     59         dzoom = dzoom * depi
     60       ELSEIF( dzoom.LT. 25. ) THEN
     61         WRITE(6,*) ' Le param. dzoomy pour fxhyp est trop petit ! L aug
     62     ,menter et relancer ! '
     63         STOP 1
     64       ELSE
     65         dzoom = dzoom * pi/180.
     66       ENDIF
     67
     68       WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)'
     69       WRITE(6,24) xzoom,grossism,tau,dzoom
    4970
    5071       DO i = 0, nmax2
    51          xtild(i) = FLOAT(i) /nmax2
    52         IF( xtild(i).EQ. 0.5 )  xtild(i) = xtild(i) + 1.e-6
    53        ENDDO
    54 
    55        DO i = 1, nmax
    56         fhyp(i) = TANH ( ( xtild(i) - 0.5*(1.- dzoom) )          /
    57      ,                 ( tau * xtild(i) * ( 0.5 -xtild(i))) )
    58        ENDDO
    59 
    60         fhyp(   0  ) = - 1.
    61         fhyp( nmax ) =   1.
     72        xtild(i) = - pi + FLOAT(i) * depi /nmax2
     73       ENDDO
     74
     75       DO i = 0, nmax2
     76
     77       fa  = tau*  ( dzoom/2.  - xtild(i) )
     78       fb  = xtild(i) *  ( pi - xtild(i) )
     79
     80
     81       IF( 200.* fb .LT. - fa )   THEN
     82         fhyp ( i) = - 1.
     83       ELSEIF( 200. * fb .LT. fa ) THEN
     84         fhyp ( i) =   1.
     85       ELSE
     86         fhyp(i) =  TANH ( fa/fb )
     87       ENDIF
     88
     89       IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
     90       IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
     91
     92       ENDDO
    6293
    6394cc  ....  Calcul  de  beta  ....
    64 c
    65        ffdx( 0 ) = 0.
    66 
    67        DO i = 1, nmax
    68         xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
    69         fxm     = TANH ( ( xmoy - 0.5 * ( 1. - dzoom ) )    /
    70      ,                 ( tau * xmoy * ( 0.5 -xmoy)) )
    71         ffdx(i) = ffdx(i-1) + fxm * ( xtild(i) - xtild(i-1) )
    72        ENDDO
    73 
    74         beta  = ( grossism * ffdx(nmax) - 0.5 ) / ( ffdx(nmax) - 0.5 )
     95c   ............................
     96
     97       ffdx = 0.
     98
     99       DO i = nmax +1,nmax2
     100
     101       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
     102       fa  = tau*  ( dzoom/2.  - xmoy )
     103       fb  = xmoy *  ( pi - xmoy )
     104
     105       IF( 200.* fb .LT. - fa )   THEN
     106         fxm = - 1.
     107       ELSEIF( 200. * fb .LT. fa ) THEN
     108         fxm =   1.
     109       ELSE
     110         fxm =  TANH ( fa/fb )
     111       ENDIF
     112
     113       IF ( xmoy.EQ. 0. )  fhyp(i) =  1.
     114       IF ( xmoy.EQ. pi )  fhyp(i) = -1.
     115       ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
     116
     117       ENDDO
     118
     119        beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )
     120
     121       IF( 2.*beta - grossism.LE. 0.)  THEN
     122        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
     123     ,tine fxhyp est mauvaise ! '
     124        WRITE(6,*)'Modifier les valeurs de  grossismx ,tau ou dzoomx ',
     125     , ' et relancer ! ***  '
     126        CALL ABORT
     127       ENDIF
    75128c
    76129c   .....  calcul  de  Xprimt   .....
    77130c
    78131       
    79        DO i = 0, nmax
     132       DO i = nmax, nmax2
    80133        Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
    81134       ENDDO
    82135c   
    83        DO i = 0, nmax
     136       DO i =  nmax+1, nmax2
    84137        Xprimt( nmax2 - i ) = Xprimt( i )
    85138       ENDDO
    … …  
    88141c   .....  Calcul  de  Xf     ........
    89142
    90         Xf(0) = 0.
    91        DO i = 1, nmax
    92         xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
    93         fxm     = TANH ( ( xmoy - 0.5 * ( 1. - dzoom ) )    /
    94      ,                 ( tau * xmoy * ( 0.5 -xmoy)) )
    95         xxpr(i)    = beta + ( grossism - beta ) * fxm
    96        ENDDO
    97 
    98        DO i = 1,nmax
    99          xxpr(nmax2-i+1) = xxpr(i)
     143       Xf(0) = - pi
     144
     145       DO i =  nmax +1, nmax2
     146
     147       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
     148       fa  = tau*  ( dzoom/2.  - xmoy )
     149       fb  = xmoy *  ( pi - xmoy )
     150
     151       IF( 200.* fb .LT. - fa )   THEN
     152         fxm = - 1.
     153       ELSEIF( 200. * fb .LT. fa ) THEN
     154         fxm =   1.
     155       ELSE
     156         fxm =  TANH ( fa/fb )
     157       ENDIF
     158
     159       IF ( xmoy.EQ. 0. )  fxm =  1.
     160       IF ( xmoy.EQ. pi )  fxm = -1.
     161       xxpr(i)    = beta + ( grossism - beta ) * fxm
     162
     163       ENDDO
     164
     165       DO i = nmax+1, nmax2
     166        xxpr(nmax2-i+1) = xxpr(i)
    100167       ENDDO
    101168
    … …  
    103170         Xf(i)   = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
    104171        ENDDO
    105         do i=1,nmax2
    106         xf(i)=xf(i)/xf(nmax2)
    107         enddo
    108 
    109 
    110        PRINT *,' XF ',xf(0),xf(nmax),xf(nmax2)
    111172
    112173
    … …  
    117178c     .....  xuv = 0.5  si  calcul  aux pts      U        ........
    118179c
     180      WRITE(6,18)
    119181c
    120182      DO 5000  ik = 1, 4
    … …  
    130192       ENDIF
    131193
     194      xo1   = 0.
    132195
    133196      DO 1500 i = 1, iim
    134197
    135       xlon2 = (  FLOAT(i) + xuv - 0.75) / FLOAT(iim)
    136 
    137       xo1   = 0.
    138       xi    = xlon2
    139 c
    140       DO 500 iter = 1,300
    141 
     198      xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim)
     199
     200      Xfi    = xlon2
     201c
    142202      DO 250 it =  nmax2,0,-1
    143       IF( xi.GE.xtild(it))  GO TO 350
     203      IF( Xfi.GE.Xf(it))  GO TO 350
    144204250   CONTINUE
    145205
    146206      it = 0
    147       xi = xtild(it)
    148207
    149208350   CONTINUE
    150        IF(it.EQ.nmax2)  THEN
    151         it       = nmax2 -1
    152         xf(it+1) = 1.
    153        ENDIF
    154 c  .................................................................
    155 c  ....  Interpolation entre  xi(it) et xi(it+1)   pour avoir X(xi) 
    156 c      .....           et   X'(xi)                             .....
    157 c  .................................................................
    158 
    159        xint   = ( Xf(it+1)-Xf(it) ) / ( xtild(it+1)-xtild(it) )      *
    160      +                    ( xi-xtild(it) )  +  Xf(it)
    161       Xprimin = ( Xprimt(it+1)-Xprimt(it) )/ ( xtild(it+1)-xtild(it) ) *
    162      +                    ( xi-xtild(it) )  +  Xprimt(it)
    163 
    164        xi = xi - (xint-xlon2)/Xprimin
    165 
    166       IF( ABS(xi-xo1).LE.epsilon) GO TO 550
    167       xo1 = xi
    168 c
     209
     210c    ......  Calcul de   Xf(xi)    ......
     211c
     212      xi  = xtild(it)
     213
     214      IF(it.EQ.nmax2)  THEN
     215       it       = nmax2 -1
     216       Xf(it+1) = pi
     217      ENDIF
     218c  .....................................................................
     219c
     220c   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
     221c   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )
     222c          et (Xf(it+1),xtild(it+1) )
     223
     224       CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
     225     ,                xtild(it),xtild(it+1),  a0, a1, a2, a3  )
     226
     227       Xf1     = Xf(it)
     228       Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
     229
     230       DO 500 iter = 1,300
     231        xi = xi - ( Xf1 - Xfi )/ Xprimin
     232
     233        IF( ABS(xi-xo1).LE.epsilon)  GO TO 550
     234         xo1      = xi
     235         xi2      = xi * xi
     236         Xf1      = a0 +  a1 * xi +     a2 * xi2  +     a3 * xi2 * xi
     237         Xprimin  =       a1      + 2.* a2 *  xi  + 3.* a3 * xi2
    169238500   CONTINUE
    170       PRINT *,' ***   PAS DE SOLUTION  ****  ',i,xlon2
    171         STOP 4
     239        WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter
     240          STOP 6
    172241550   CONTINUE
    173242
    174        xxprim(i)   = depi/( FLOAT(iim) * Xprimin)
    175        xvrai(i)    = depi *  (xi - 0.5) + xzoom
    176 
     243       xxprim(i) = depi/ ( FLOAT(iim) * Xprimin )
     244       xvrai(i)  =  xi + xzoom
    177245
    1782461500   CONTINUE
    179247
    180248       DO i = 1 , iim
    181         xlon  (i)   = xvrai(i)
     249        xlon(i)     = xvrai(i)
    182250        xprimm(i)   = xxprim(i)
    183 cc        xxlon(i) = xlon(i)*180./pi
    184251       ENDDO
    185252       
    186 cc      PRINT *,' XLON avant '
    187 cc      PRINT 68,(xxlon(i),i=1,iim)
    188        
    189253       DO i = 1, iim -1
    190254        IF( xvrai(i+1). LT. xvrai(i) )  THEN
    191          PRINT *,' PBS.  avec  rlonu(',i+1,' plus petit que rlonu(',i,
    192      ,   ')'
    193          STOP
     255         WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
     256     ,  ')'
     257        STOP 7
    194258        ENDIF
    195259       ENDDO
    … …  
    205269       ENDDO
    206270
    207        PRINT *,' LONGITUDES  min max ',champmin,champmax
    208 
    209271      IF(champmin .GE. - pi .AND. champmax.LE. pi )  THEN
    210272                GO TO 1600
    211273      ELSE
    212         PRINT 18
    213         PRINT *,'Reorganisation des longitudes pour avoir entre - pi ',
    214      ,  ' et  pi '
     274       WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
     275     ,  ' et pi '
    215276c
    216277        IF( xzoom.LE.0.)  THEN
    … …  
    219280           IF( xvrai(i).GE. - pi )  GO TO 80
    220281          ENDDO
    221             PRINT *, ' PBS. 1 '
    222             STOP
     282            WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
     283            STOP 8
    223284 80       CONTINUE
    224285          is2 = i
    … …  
    240301            IF( xvrai(i).LE. pi )  GO TO 90
    241302           ENDDO
    242              PRINT *,' PBS.  2 '
    243               STOP
     303             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
     304              STOP 9
    244305 90        CONTINUE
    245306            is2 = i
    246307          ENDIF
    247 cc           PRINT *,' IS2 ',is2
    248308           idif = iim -is2
    249309           DO ii = 1, is2
    … …  
    270330         ENDDO
    271331
    272 
    273332         IF( ik.EQ.1 )  THEN
    274          PRINT *, ' XLON aux pts. V-0.25   apres ( en  deg. ) '
    275          PRINT 18
    276          PRINT 68,xvrai
    277          PRINT *,' XPRIM '
    278          PRINT 68, xprimm
     333c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
     334c          WRITE(6,18)
     335c          WRITE(6,68) xvrai
     336ccc          WRITE(6,*) ' XPRIM k ',ik
     337ccc          WRITE(6,566)  xprimm
     338
    279339           DO i = 1,iim + 1
    280340             rlonm025(i) = xlon( i )
    281341            xprimm025(i) = xprimm(i)
    282342           ENDDO
     343
    283344         ELSE IF( ik.EQ.2 )  THEN
    284          PRINT 18
    285          PRINT *, ' XLON aux pts. V   apres ( en  deg. ) '
    286          PRINT 68,xvrai
    287          PRINT *,' XPRIM '
    288          PRINT 68, xprimm
     345c          WRITE(6,18)
     346c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
     347c          WRITE(6,68) xvrai
     348ccc          WRITE(6,*) ' XPRIM k ',ik
     349ccc          WRITE(6,566)  xprimm
     350
    289351           DO i = 1,iim + 1
    290352             rlonv(i) = xlon( i )
    291353            xprimv(i) = xprimm(i)
    292354           ENDDO
    293          ELSE IF( ik.EQ.3 )  THEN
    294          PRINT 18
    295          PRINT *, ' XLON aux pts. U   apres ( en  deg. ) '
    296          PRINT 68,xvrai
    297          PRINT *,' XPRIM '
    298          PRINT 68, xprimm
     355
     356         ELSE IF( ik.EQ.3)   THEN
     357c          WRITE(6,18)
     358c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
     359c          WRITE(6,68) xvrai
     360ccc          WRITE(6,*) ' XPRIM ik ',ik
     361ccc          WRITE(6,566)  xprimm
     362
    299363           DO i = 1,iim + 1
    300364             rlonu(i) = xlon( i )
    301365            xprimu(i) = xprimm(i)
    302366           ENDDO
     367
    303368         ELSE IF( ik.EQ.4 )  THEN
    304          PRINT 18
    305          PRINT *, ' XLON aux pts. V+0.25   apres ( en  deg. ) '
    306          PRINT 68,xvrai
    307          PRINT *,' XPRIM '
    308          PRINT 68, xprimm
     369c          WRITE(6,18)
     370c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
     371c          WRITE(6,68) xvrai
     372ccc          WRITE(6,*) ' XPRIM ik ',ik
     373ccc          WRITE(6,566)  xprimm
     374
    309375           DO i = 1,iim + 1
    310376             rlonp025(i) = xlon( i )
    311377            xprimp025(i) = xprimm(i)
    312378           ENDDO
     379
    313380         ENDIF
    314381
    3153825000    CONTINUE
    316383c
    317 c       ...........  fin  de la boucle  do 5000      ............
    318 
    319 c
    320         DO i = 1, iim + 1
    321          dlon1(i) = rlonm025(i) - rlonv(i)
    322          dlon2(i) = rlonm025(i) - rlonp025(i)
    323          dlon3(i) = rlonm025(i) - rlonu(i)
     384       WRITE(6,18)
     385c
     386c    ...........  fin  de la boucle  do 5000      ............
     387
     388        DO i = 1, iim
     389         xlon(i) = rlonv(i+1) - rlonv(i)
    324390        ENDDO
    325 
    326         DO i = 1, iim + 1
    327          rlonm025(i) = rlonm025(i) + dlon1(i)
     391        champmin =  1.e12
     392        champmax = -1.e12
     393        DO i = 1, iim
     394         champmin = MIN( champmin, xlon(i) )
     395         champmax = MAX( champmax, xlon(i) )
    328396        ENDDO
    329 
    330         DO i = 1, iim + 1
    331          rlonv(i)    = rlonm025(i) - dlon1(i)
    332          rlonp025(i) = rlonm025(i) - dlon2(i)
    333          rlonu(i)    = rlonm025(i) - dlon3(i)
    334         ENDDO
    335 
    336         DO i = 1, iim
    337          xprimu   (i) = rlonu(i+1) - rlonu(i)
    338          xprimv   (i) = rlonv(i+1) - rlonv(i)
    339          xprimm025(i) = rlonm025(i+1) - rlonm025(i)
    340          xprimp025(i) = rlonp025(i+1) - rlonp025(i)
    341         ENDDO
    342          xprimu   (iip1) = xprimu   (1)
    343          xprimv   (iip1) = xprimv   (1)
    344          xprimm025(iip1) = xprimm025(1)
    345          xprimp025(iip1) = xprimp025(1)
    346 
    347 
    348 18       FORMAT(/)
    349 68       FORMAT(1x,7f9.2)
    350 
    351 
    352          RETURN
    353          END
     397         champmin = champmin * 180./pi
     398         champmax = champmax * 180./pi
     399
     40018     FORMAT(/)
     40124     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
     40268     FORMAT(1x,7f9.2)
     403
     404       RETURN
     405       END

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxyhyper.F

    r2 r232  
    1        SUBROUTINE fxyhyper ( yzoom, grossy, dzoomy,tauy,deltay  ,
    2      ,                      xzoom, grossx, dzoomx,taux,
     1c
     2c $Header$
     3c
     4       SUBROUTINE fxyhyper ( yzoom, grossy, dzoomy,tauy  ,   
     5     ,                       xzoom, grossx, dzoomx,taux  ,
    36     , rlatu,yprimu,rlatv,yprimv,rlatu1,  yprimu1,  rlatu2,  yprimu2  ,
    47     , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
    … …  
    1922c     a) le grossissement du zoom  :  grossy  ( en y ) et grossx ( en x )
    2023c     b) l' extension     du zoom  :  dzoomy  ( en y ) et dzoomx ( en x )
    21 c     c) la raideur       du zoom  :   tau   , ici  =  1.
     24c     c) la raideur de la transition du zoom  :   taux et tauy   
    2225c
    23 c     N.B :   le produit ( grossissement * extension )  doit etre <  1.
    24 c    *******
    25 c     En plus , il y a un autre parametre , moins important mais quand
    26 c     meme utile , c'est  deltay , deplacement de la zone du zoom en
    27 c     latitude   :  Si  on deplace  de 10. degres vers le nord  ( deltay
    28 c      = 10.  dans inigeom ) .
    29 c
     26c  N.B : Il vaut mieux avoir   :   grossx * dzoomx <  pi    ( radians )
     27c ******
     28c                  et              grossy * dzoomy <  pi/2  ( radians )
    3029c
    3130#include "dimensions.h"
    … …  
    4039       REAL rlonu(iip1),xprimu(iip1),rlonv(iip1),xprimv(iip1),
    4140     , rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),xprimp025(iip1)
    42        REAL deltay
     41       REAL*8  dxmin, dxmax , dymin, dymax
    4342
    4443c   ....   var. locales   .....
    … …  
    4746c
    4847
    49        CALL fyhyp ( yzoom, grossy, dzoomy,tauy, deltay ,
    50      ,  rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 )
     48       CALL fyhyp ( yzoom, grossy, dzoomy,tauy  ,
     49     ,  rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
     50     ,  dymin,dymax                                               )
     51
     52       CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv,
     53     , xprimv,rlonu,xprimu,rlonp025,xprimp025 , dxmin,dxmax         )
    5154
    5255
    53        CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv,
    54      , xprimv,rlonu,xprimu,rlonp025,xprimp025 )
    55 
    56 
    57         DO i = 1, iim
     56        DO i = 1, iip1
    5857          IF(rlonp025(i).LT.rlonv(i))  THEN
    59            PRINT *,' Attention !  rlonp025 < rlonv',i
     58           WRITE(6,*) ' Attention !  rlonp025 < rlonv',i
    6059            STOP
    6160          ENDIF
    6261
    6362          IF(rlonv(i).LT.rlonm025(i))  THEN
    64            PRINT *,' Attention !  rlonm025 > rlonv',i
     63           WRITE(6,*) ' Attention !  rlonm025 > rlonv',i
    6564            STOP
    6665          ENDIF
    6766
    6867          IF(rlonp025(i).GT.rlonu(i))  THEN
    69            PRINT *,' Attention !  rlonp025 > rlonu',i
     68           WRITE(6,*) ' Attention !  rlonp025 > rlonu',i
    7069            STOP
    7170          ENDIF
    7271        ENDDO
    7372
    74         PRINT *,'  *** TEST DE COHERENCE  OK    POUR   FX **** '
     73        WRITE(6,*) '  *** TEST DE COHERENCE  OK    POUR   FX **** '
    7574
    7675c
    … …  
    7877c
    7978       IF(rlatu1(j).LE.rlatu2(j))   THEN
    80          PRINT *,' Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j
     79         WRITE(6,*)'Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j
    8180         STOP 13
    8281       ENDIF
    8382c
    8483       IF(rlatu2(j).LE.rlatu(j+1))  THEN
    85         PRINT *,' Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j
     84        WRITE(6,*)'Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j
    8685        STOP 14
    8786       ENDIF
    8887c
    8988       IF(rlatu(j).LE.rlatu1(j))    THEN
    90         PRINT *,' Attention ! rlatu < rlatu1 ',rlatu(j),rlatu1(j),j
     89        WRITE(6,*)' Attention ! rlatu < rlatu1 ',rlatu(j),rlatu1(j),j
    9190        STOP 15
    9291       ENDIF
    9392c
    9493       IF(rlatv(j).LE.rlatu2(j))    THEN
    95         PRINT *,' Attention ! rlatv < rlatu2 ',rlatv(j),rlatu2(j),j
     94        WRITE(6,*)' Attention ! rlatv < rlatu2 ',rlatv(j),rlatu2(j),j
    9695        STOP 16
    9796       ENDIF
    9897c
    9998       IF(rlatv(j).ge.rlatu1(j))    THEN
    100         PRINT *,' Attention ! rlatv > rlatu1 ',rlatv(j),rlatu1(j),j
     99        WRITE(6,*)' Attention ! rlatv > rlatu1 ',rlatv(j),rlatu1(j),j
    101100        STOP 17
    102101       ENDIF
    103102c
    104103       IF(rlatv(j).ge.rlatu(j))     THEN
    105         PRINT *,'Attention ! rlatv > rlatu ',rlatv(j),rlatu(j),j
     104        WRITE(6,*) ' Attention ! rlatv > rlatu ',rlatv(j),rlatu(j),j
    106105        STOP 18
    107106       ENDIF
    … …  
    109108       ENDDO
    110109c
    111        PRINT *,'  *** TEST DE COHERENCE  OK    POUR   FY **** '
    112 
     110       WRITE(6,*) '  *** TEST DE COHERENCE  OK    POUR   FY **** '
    113111c
     112        WRITE(6,18)
     113        WRITE(6,*) '  Latitudes  '
     114        WRITE(6,*) ' *********** '
     115        WRITE(6,18)
     116        WRITE(6,3)  dymin, dymax
     117        WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par
     118     ,ametres  grossism , tau , dzoom pour Y et repasser ! '
     119c
     120        WRITE(6,18)
     121        WRITE(6,*) '  Longitudes  '
     122        WRITE(6,*) ' ************ '
     123        WRITE(6,18)
     124        WRITE(6,3)  dxmin, dxmax
     125        WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par
     126     ,ametres  grossism , tau , dzoom pour Y et repasser ! '
     127        WRITE(6,18)
     128c
     1293      Format(1x, ' Au centre du zoom , la longueur de la maille est',
     130     ,  ' d environ ',f8.2 ,' degres  ',
     131     , ' alors que la maille en dehors de la zone du zoom est d environ
     132     , ', f8.2,' degres ' )
     13318      FORMAT(/)
    114134
    115135       RETURN

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/fyhyp.F

    r2 r232  
    1        SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau,deltay , 
    2      ,  rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 )
    3 
     1c
     2c $Header$
     3c
     4       SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau  , 
     5     ,  rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
     6     ,  champmin,champmax                                            )
     7
     8cc    ...  Version du 01/04/2001 ....
    49
    510       IMPLICIT NONE
    … …  
    1318c
    1419c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc)
    15 c     dzoom  etant  la distance totale de la zone du zoom
    16 c     tau  la transition ,   normalement  = 1  .
    17 
    18 c  N.B :  on doit avoir :  grossism * dzoom  <   1 .
    19 c         **************
    20 c
    21 c    Pour Indoex , on a pris :
    22 c         *******
    23 c    grossism = 2.5 , dzoom = 7/24  en x et  y  , pour iim = 128 et jjm=64
    24 c    yzoomdeg = 0.  , tau = 1.  et delaty = 10.
     20c     dzoom  etant  la distance totale de la zone du zoom ( en radians )
     21c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom   
     22c
     23c
     24c N.B : Il vaut mieux avoir : grossism * dzoom  <  pi/2  (radians) ,en lati.
     25c      ********************************************************************
    2526c
    2627c
    … …  
    2930
    3031       INTEGER      nmax , nmax2
    31        PARAMETER (  nmax = 50000, nmax2 = 2*nmax )
     32       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
    3233c
    3334c
    3435c     .......  arguments  d'entree    .......
    3536c
    36        REAL yzoomdeg, grossism,dzoom,tau , deltay
     37       REAL yzoomdeg, grossism,dzoom,tau
     38c         ( rentres  par  run.def )
    3739
    3840c     .......  arguments  de sortie   .......
    … …  
    4244
    4345c
    44 c    ..... Champs  locaux    .....
     46c     .....     champs  locaux    .....
    4547c
    4648     
    47        REAl ylat(jjp1), yprim(jjp1)
    48        REAL yuv
    49        REAL ytild(0:nmax2)
    50        REAL fhyp(0:nmax),ffdx(0:nmax),beta,Ytprim(0:nmax2)
    51        SAVE Ytprim, ytild,Yf
    52        REAL Yf(0:nmax2),yypr(0:nmax2)
    53        REAL yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
    54        REAL pi,depi,pis2,epsilon,yzoom
    55        REAL yo1,yi,ylon2,fxm,ymoy,yint,Yprimin
    56        REAL ypn,deply,y00
     49       REAL*8 ylat(jjp1), yprim(jjp1)
     50       REAL*8 yuv
     51       REAL*8 yt(0:nmax2)
     52       REAL*8 fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2)
     53       SAVE Ytprim, yt,Yf
     54       REAL*8 Yf(0:nmax2),yypr(0:nmax2)
     55       REAL*8 yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
     56       REAL*8 pi,depi,pis2,epsilon,y0,pisjm
     57       REAL*8 yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax
     58       REAL*8 yfi,Yf1,ffdy
     59       REAL*8 ypn,deply,y00
    5760       SAVE y00, deply
    5861
    … …  
    6063       INTEGER jpn,jjpn
    6164       SAVE jpn
    62 
     65       REAL*8 a0,a1,a2,a3,yi2,heavyy0,heavyy0m
     66       REAL*8 fa(0:nmax2),fb(0:nmax2)
     67       REAL y0min,y0max
     68
     69       REAL*8     heavyside
     70       EXTERNAL   heavyside
    6371
    6472       pi       = 2. * ASIN(1.)
    6573       depi     = 2. * pi
    6674       pis2     = pi/2.
    67        epsilon  = 1.e-6
    68        yzoom    = yzoomdeg * pi/180.
    69 
    70 
     75       pisjm    = pi/ FLOAT(jjm)
     76       epsilon  = 1.e-3
     77       y0       =  yzoomdeg * pi/180.
     78
     79       IF( dzoom.LT.1.)  THEN
     80         dzoom = dzoom * pi
     81       ELSEIF( dzoom.LT. 12. ) THEN
     82         WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug
     83     ,menter et relancer ! '
     84         STOP 1
     85       ELSE
     86         dzoom = dzoom * pi/180.
     87       ENDIF
     88
     89       WRITE(6,18)
     90       WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)'
     91       WRITE(6,24) y0,grossism,tau,dzoom
    7192
    7293       DO i = 0, nmax2
    73         ytild(i) = FLOAT(i) /nmax2
    74        IF( ytild(i).EQ.0.5 )  ytild(i) = ytild(i) + 1.e-6
    75        ENDDO
    76 
    77        DO i = 1, nmax
    78         fhyp(i) = TANH ( ( ytild(i) - 0.5*(1.- dzoom) )          /
    79      ,                 ( tau * ytild(i) * ( 0.5 -ytild(i))) )
    80        ENDDO
    81 
    82        fhyp(   0  ) = - 1.
    83        fhyp( nmax ) =   1.
     94        yt(i) = - pis2  + FLOAT(i)* pi /nmax2
     95       ENDDO
     96
     97       heavyy0m = heavyside( -y0 )
     98       heavyy0  = heavyside(  y0 )
     99       y0min    = 2.*y0*heavyy0m - pis2
     100       y0max    = 2.*y0*heavyy0  + pis2
     101
     102       DO i = 0, nmax2
     103        IF( yt(i).LT.y0 )  THEN
     104         fa (i) = tau*  (yt(i)-y0+dzoom/2. )
     105         fb(i) =   (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) )
     106        ELSEIF ( yt(i).GT.y0 )  THEN
     107         fa(i) =   tau *(y0-yt(i)+dzoom/2. )
     108         fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 )
     109       ENDIF
     110       
     111       IF( 200.* fb(i) .LT. - fa(i) )   THEN
     112         fhyp ( i) = - 1.
     113       ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
     114         fhyp ( i) =   1.
     115       ELSE 
     116         fhyp(i) =  TANH ( fa(i)/fb(i) )
     117       ENDIF
     118
     119       IF( yt(i).EQ.y0 )  fhyp(i) = 1.
     120       IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1.
     121
     122       ENDDO
    84123
    85124cc  ....  Calcul  de  beta  ....
    86125c
    87        ffdx( 0 ) = 0.
    88 
    89        DO i = 1, nmax
    90         ymoy    = 0.5 * ( ytild(i-1) + ytild( i ) )
    91         fxm     = TANH ( ( ymoy - 0.5 * ( 1. - dzoom ) )    /
    92      ,                 ( tau * ymoy * ( 0.5 -ymoy)) )
    93         ffdx(i) = ffdx(i-1) + fxm * ( ytild(i) - ytild(i-1) )
    94        ENDDO
    95 
    96         beta  = ( grossism * ffdx(nmax) - 0.5 ) / ( ffdx(nmax) - 0.5 )
     126       ffdy   = 0.
     127
     128       DO i = 1, nmax2
     129        ymoy    = 0.5 * ( yt(i-1) + yt( i ) )
     130        IF( ymoy.LT.y0 )  THEN
     131         fa(i)= tau * ( ymoy-y0+dzoom/2.)
     132         fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy )
     133        ELSEIF ( ymoy.GT.y0 )  THEN
     134         fa(i)= tau * ( y0-ymoy+dzoom/2. )
     135         fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 )
     136        ENDIF
     137
     138        IF( 200.* fb(i) .LT. - fa(i) )    THEN
     139         fxm ( i) = - 1.
     140        ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
     141         fxm ( i) =   1.
     142        ELSE
     143         fxm(i) =  TANH ( fa(i)/fb(i) )
     144        ENDIF
     145         IF( ymoy.EQ.y0 )  fxm(i) = 1.
     146         IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1.
     147         ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) )
     148
     149        ENDDO
     150
     151        beta  = ( grossism * ffdy - pi ) / ( ffdy - pi )
     152
     153       IF( 2.*beta - grossism.LE. 0.)  THEN
     154
     155        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
     156     ,tine fyhyp est mauvaise ! '
     157        WRITE(6,*)'Modifier les valeurs de  grossismy ,tauy ou dzoomy',
     158     , ' et relancer ! ***  '
     159        CALL ABORT
     160
     161       ENDIF
    97162c
    98163c   .....  calcul  de  Ytprim   .....
    99164c
    100165       
    101        DO i = 0, nmax
     166       DO i = 0, nmax2
    102167        Ytprim(i) = beta  + ( grossism - beta ) * fhyp(i)
    103168       ENDDO
    104 c   
    105        DO i = 0, nmax
    106         Ytprim( nmax2 - i ) = Ytprim( i )
    107        ENDDO
    108 c
    109169
    110170c   .....  Calcul  de  Yf     ........
    111171
    112         Yf(0) = 0.
    113        DO i = 1, nmax
    114         ymoy    = 0.5 * ( ytild(i-1) + ytild( i ) )
    115         fxm     = TANH ( ( ymoy - 0.5 * ( 1. - dzoom ) )    /
    116      ,                 ( tau * ymoy * ( 0.5 -ymoy)) )
    117         yypr(i)    = beta + ( grossism - beta ) * fxm
    118        ENDDO
    119 
    120        DO i = 1,nmax
    121          yypr(nmax2-i+1) = yypr(i)
    122        ENDDO
    123 
    124         DO i=1,nmax2
    125          Yf(i)   = Yf(i-1) + yypr(i) * ( ytild(i) - ytild(i-1) )
    126         ENDDO
    127 c
    128 c
     172       Yf(0) = - pis2
     173       DO i = 1, nmax2
     174        yypr(i)    = beta + ( grossism - beta ) * fxm(i)
     175       ENDDO
     176
     177       DO i=1,nmax2
     178        Yf(i)   = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) )
     179       ENDDO
    129180
    130181c    ****************************************************************
    … …  
    133184c   .....   yuv  = 0.5  si calcul des latitudes  aux pts.  V  .....
    134185c
     186      WRITE(6,18)
    135187c
    136188      DO 5000  ik = 1,4
    … …  
    149201         jlat = jjm
    150202       ENDIF
    151        
    152 c
     203c
     204       yo1   = 0.
    153205       DO 1500 j =  1,jlat
    154 
    155         ylon2 =  ( FLOAT(j)  + yuv  -1.) / FLOAT(jjm)
    156 
    157206        yo1   = 0.
    158         yi    = ylon2
    159 
    160 c
    161        DO 500 iter = 1,300
    162 
     207        ylon2 =  - pis2 + pisjm * ( FLOAT(j)  + yuv  -1.) 
     208        yfi    = ylon2
     209c
    163210       DO 250 it =  nmax2,0,-1
    164         IF( yi.GE.ytild(it))  GO TO 350
     211        IF( yfi.GE.Yf(it))  GO TO 350
    165212250    CONTINUE
    166 
    167213       it = 0
    168        yi = ytild(it)
    169 
    170214350    CONTINUE
    171215
     216       yi = yt(it)
    172217       IF(it.EQ.nmax2)  THEN
    173218        it       = nmax2 -1
    174         Yf(it+1) = 1.
     219        Yf(it+1) = pis2
    175220       ENDIF
    176221c  .................................................................
    … …  
    179224c  .................................................................
    180225
    181        yint   = ( Yf(it+1)-Yf(it) ) / ( ytild(it+1)-ytild(it) )      *
    182      +                    ( yi-ytild(it) )  +  Yf(it)
    183       Yprimin = ( Ytprim(it+1)-Ytprim(it) )/ ( ytild(it+1)-ytild(it) ) *
    184      +                    ( yi-ytild(it) )  +  Ytprim(it)
    185        yi = yi - (yint-ylon2)/Yprimin
    186 
    187       IF( ABS(yi-yo1).LE.epsilon) GO TO 550
    188       yo1 = yi
    189 c
     226       CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1),   
     227     ,                  yt(it),yt(it+1) ,   a0,a1,a2,a3   )     
     228
     229       Yf1     = Yf(it)
     230       Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi
     231
     232       DO 500 iter = 1,300
     233         yi = yi - ( Yf1 - yfi )/ Yprimin
     234
     235        IF( ABS(yi-yo1).LE.epsilon)  GO TO 550
     236         yo1      = yi
     237         yi2      = yi * yi
     238         Yf1      = a0 +  a1 * yi +     a2 * yi2  +     a3 * yi2 * yi
     239         Yprimin  =       a1      + 2.* a2 *  yi  + 3.* a3 * yi2
    190240500   CONTINUE
    191       PRINT *,' ***   PAS DE SOLUTION  **** ',j,ylon2,iter
    192         STOP 4
     241        WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter
     242         STOP 2
    193243550   CONTINUE
    194 
    195        yprim(j)    = pi /( FLOAT(jjm) *  Yprimin)
    196        yvrai(j)    = pi *  (yi - 0.5) + yzoom
     244c
     245       Yprimin   = a1  + 2.* a2 *  yi   + 3.* a3 * yi* yi
     246       yprim(j)  = pi / ( jjm * Yprimin )
     247       yvrai(j)  = yi
    197248
    1982491500    CONTINUE
    199 
    200 cc          print *,' LAT avant reorgan '
    201 cc         print 68,(yyvrai(j),j=1,jlat)
    202250
    203251       DO j = 1, jlat -1
    204252        IF( yvrai(j+1). LT. yvrai(j) )  THEN
    205          PRINT *,' PBS.  avec  rlat(',j+1,' plus petit que rlat(',j,
    206      ,   ')'
    207          STOP
    208         ENDIF
    209        ENDDO
    210 
    211         PRINT 18
    212         PRINT *,'Reorganisation des latitudes pour avoir entre - pi/2 ',
    213      ,  ' et  pi/2 '
    214 c
    215            
     253         WRITE(6,*) ' PBS. avec  rlat(',j+1,') plus petit que rlat(',j,
     254     ,  ')'
     255         STOP 3
     256        ENDIF
     257       ENDDO
     258
     259       WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2'
     260     , ,' et  pi/2 '
     261c
    216262        IF( ik.EQ.1 )   THEN
    217            ypn = pis2 - deltay * pi/180.
     263           ypn = pis2
    218264          DO j = jlat,1,-1
    219265           IF( yvrai(j).LE. ypn ) GO TO 1502
    … …  
    239285         ENDDO
    240286
    241 
    242 
    243287c      ***********   Fin de la reorganisation     *************
    244288c
    245289 1600   CONTINUE
    246 
    247 
    248290
    249291       DO j = 1, jlat
    … …  
    256298        ENDDO
    257299
    258 
    259300        IF( ik.EQ.1 )  THEN
    260         PRINT 18
    261         PRINT *, ' YLAT  en U   apres ( en  deg. ) '
    262         PRINT 68,(yvrai(j),j=1,jlat)
    263         PRINT *,' YPRIM '
    264         PRINT 68,( yprim(j),j=1,jlat)
     301c         WRITE(6,18)
     302c         WRITE(6,*)  ' YLAT  en U   apres ( en  deg. ) '
     303c         WRITE(6,68) (yvrai(j),j=1,jlat)
     304cc         WRITE(6,*) ' YPRIM '
     305cc         WRITE(6,445) ( yprim(j),j=1,jlat)
     306
    265307          DO j = 1, jlat
    266308            rrlatu(j) =  ylat( j )
    267309           yyprimu(j) = yprim( j )
    268310          ENDDO
    269 c
     311
    270312        ELSE IF ( ik.EQ. 2 )  THEN
    271         PRINT 18
    272         PRINT *, ' YLAT   en V  apres ( en  deg. ) '
    273         PRINT 68,(yvrai(j),j=1,jlat)
    274         PRINT *,' YPRIM '
    275         PRINT 68,( yprim(j),j=1,jlat)
     313c         WRITE(6,18)
     314c         WRITE(6,*) ' YLAT   en V  apres ( en  deg. ) '
     315c         WRITE(6,68) (yvrai(j),j=1,jlat)
     316cc         WRITE(6,*)' YPRIM '
     317cc         WRITE(6,445) ( yprim(j),j=1,jlat)
     318
    276319          DO j = 1, jlat
    277320            rrlatv(j) =  ylat( j )
    278321           yyprimv(j) = yprim( j )
    279322          ENDDO
    280 c
     323
    281324        ELSE IF ( ik.EQ. 3 )  THEN
    282         PRINT 18
    283         PRINT *, ' YLAT  en U + 0.75  apres ( en  deg. ) '
    284         PRINT 68,(yvrai(j),j=1,jlat)
    285         PRINT *,' YPRIM '
    286         PRINT 68,( yprim(j),j=1,jlat)
     325c         WRITE(6,18)
     326c         WRITE(6,*)  ' YLAT  en U + 0.75  apres ( en  deg. ) '
     327c         WRITE(6,68) (yvrai(j),j=1,jlat)
     328cc         WRITE(6,*) ' YPRIM '
     329cc         WRITE(6,445) ( yprim(j),j=1,jlat)
     330
    287331          DO j = 1, jlat
    288332            rlatu2(j) =  ylat( j )
    … …  
    291335
    292336        ELSE IF ( ik.EQ. 4 )  THEN
    293         PRINT 18
    294         PRINT *, ' YLAT en U + 0.25  apres ( en  deg. ) '
    295         PRINT 68,(yvrai(j),j=1,jlat)
    296         PRINT *,' YPRIM '
    297         PRINT 68,( yprim(j),j=1,jlat)
     337c         WRITE(6,18)
     338c         WRITE(6,*)  ' YLAT en U + 0.25  apres ( en  deg. ) '
     339c         WRITE(6,68)(yvrai(j),j=1,jlat)
     340cc         WRITE(6,*) ' YPRIM '
     341cc         WRITE(6,68) ( yprim(j),j=1,jlat)
     342
    298343          DO j = 1, jlat
    299344            rlatu1(j) =  ylat( j )
    300345           yprimu1(j) = yprim( j )
    301346          ENDDO
     347
    302348        ENDIF
    303349
    3043505000   CONTINUE
    305351c
     352        WRITE(6,18)
     353c
    306354c  .....     fin de la boucle  do 5000 .....
    307355
     356        DO j = 1, jjm
     357         ylat(j) = rrlatu(j) - rrlatu(j+1)
     358        ENDDO
     359        champmin =  1.e12
     360        champmax = -1.e12
     361        DO j = 1, jjm
     362         champmin = MIN( champmin, ylat(j) )
     363         champmax = MAX( champmax, ylat(j) )
     364        ENDDO
     365         champmin = champmin * 180./pi
     366         champmax = champmax * 180./pi
     367
     36824     FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
    30836918      FORMAT(/)
    30937068      FORMAT(1x,7f9.2)

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F

    r228 r232  
    2222c
    2323c=======================================================================
    24 
    25 c   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv
    26 c        et possibilite d'appeler une fonction f(y)  a derivee tangente
    27 c        hyperbolique a la  place de la fonction a derivee sinusoidale.         
    28 
    29 c   ...  Possibilite de choisir le shema de Van-leer pour l'advection de
    30 c         q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
     24c
     25c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
     26c      et possibilite d'appeler une fonction f(y)  a derivee tangente
     27c      hyperbolique a la  place de la fonction a derivee sinusoidale.
     28
     29c  ... Possibilite de choisir le shema de Van-leer pour l'advection de
     30c        q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
     31c
     32c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
    3133c
    3234c-----------------------------------------------------------------------
    … …  
    135137      LOGICAL true_calendar
    136138      PARAMETER (true_calendar = .false.)
     139C Run nudge
     140      LOGICAL ok_nudge
     141      PARAMETER (ok_nudge = .false.)
    137142
    138143c-----------------------------------------------------------------------
    … …  
    160165c     iadv = 2    shema   amont
    161166c     iadv = 3    shema  Van-leer
     167c     iadv = 4    schema  Van-leer + humidite specifique
     168c                        Modif F.Codron
    162169c
    163170c     dans le tableau q(ij,l,iq) , iq = 1  pour l'eau vapeur
    … …  
    171178      ENDDO
    172179c
    173       DO  iq = 1, nqmx
     180       IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
     181     * ,' pour la vapeur d''eau'
     182       IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema  amont',' pour la'
     183     * ,' vapeur d''eau '
     184       IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
     185     * ,'la vapeur d''eau'
     186       IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema  Van-Leer + humidite'
     187     * ,' specifique pour la vapeur d''eau'
     188c
     189       IF( iadv(1).LE.0.OR.iadv(1).GT.4 )   THEN
     190        PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
     191     * ,  car  iadv(1) = ', iadv(1)
     192         CALL ABORT
     193       ENDIF
     194
     195      DO  iq = 2, nqmx
    174196       IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
    175197     * ,' pour le traceur no ', iq
    … …  
    178200       IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
    179201     * ,'le traceur no ', iq
    180        IF( iadv(iq).LE.0.OR.iadv(iq).GT.3 )   THEN
     202
     203       IF( iadv(iq).EQ.4 )  THEN
     204         PRINT *,' Le shema  Van-Leer + humidite specifique ',
     205     * ' est  uniquement pour la vapeur d eau .'
     206         PRINT *,' Corriger iadv( ',iq, ')  et repasser ! '
     207         CALL ABORT
     208       ENDIF
     209
     210       IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 )   THEN
    181211       PRINT *,' Erreur dans le  choix de  iadv . Corriger et repasser
    182212     * . '
    183          STOP
     213         CALL ABORT
    184214       ENDIF
    185215      ENDDO
    … …  
    188218      numvanle = nqmx + 1
    189219      DO  iq = 1, nqmx
    190         IF( iadv(iq).EQ.3.AND.first ) THEN
     220        IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
    191221          numvanle = iq
    192222          first    = .FALSE.
    … …  
    195225c
    196226      DO  iq = 1, nqmx
    197         IF( iadv(iq).NE.3.AND.iq.GT.numvanle )   THEN
     227      IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle )  THEN
    198228          PRINT *,' Il y a discontinuite dans le choix du shema de ',
    199229     *    'Van-leer pour les traceurs . Corriger et repasser . '
    200            STOP
    201         ENDIF
    202         IF( iadv(iq).LT.1.OR.iadv(iq).GT.3 )    THEN
     230           CALL ABORT
     231      ENDIF
     232        IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 )    THEN
    203233          PRINT *,' Le choix de  iadv  est errone pour le traceur ',
    204234     *    iq
    … …  
    297327
    298328   1  CONTINUE
    299 c     call nudge(itau,ucov,vcov,teta,masse,ps)
     329
     330      if (ok_nudge) then
     331        call nudge(itau,ucov,vcov,teta,masse,ps)
     332      endif
    300333c
    301334c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
    … …  
    379412c
    380413            CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
    381      *                             p,  masse, dq                  )
     414     *                      p, masse, dq,  iadv(1), teta, pk     )
     415c
     416c                   ...  Modif  F.Codron  ....
    382417c
    383418         ENDIF
    … …  
    389424
    390425            CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
    391      .   time_step, itau,fluxid, fluxvid,fluxdid )
     426     .   time_step, itau)
     427
    392428
    393429         ENDIF

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/inigeom.F

    r2 r232  
     1c
     2c $Header$
     3c
    14      SUBROUTINE inigeom
    25c
    36c     Auteur :  P. Le Van
    4 c    .....................
    5 c
    6 c   ............      Version  du 20/12/98     ........................
     7c
     8c   ............      Version  du 01/04/2001     ........................
    79c
    810c  Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4  aux memes en-
    911c     endroits que les aires aireij1,..aireij4 .
     12
    1013c  Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
    1114c
    1215c
    1316      IMPLICIT NONE
    14 c
    15       REAL        deltay, tauy, taux
    16       PARAMETER ( deltay =  0. , tauy = 1. , taux = 1. )
    17 c
    18 c     deltay  est ( en degres ) le deplacement eventuel en Y du zoom
    19 cc    taux  et  tauy  sont  les raideurs   du  zoom
    20 c
    2117c
    2218#include "dimensions.h"
    … …  
    5046      SAVE rlonm025,xprimm025,rlonp025,xprimp025
    5147
    52 
    53 
    54 c----------------------------------------------------------------------
    5548      REAL      SSUM
    5649      EXTERNAL  SSUM
    57 c
    5850c
    5951c
    … …  
    6254c   -    calcul des coeff. ( cu, cv , 1./cu**2,  1./cv**2  )         -
    6355c   -                                                                -
    64 c   ------------------------------------------------------------------
    6556c   ------------------------------------------------------------------
    6657c
    … …  
    168159c
    169160c
    170       PRINT 3
     161      WRITE(6,3)
    171162 3    FORMAT( // 10x,' ....  INIGEOM  date du 01/06/98   ..... ',
    172163     * //5x,'   Calcul des elongations cu et cv  comme sommes des 4 ' /
    … …  
    191182      ENDIF
    192183
    193       PRINT *,' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
     184      WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
    194185     *  nitergdiv,nitergrot,niterh
    195186c
    196187      pi    = 2.* ASIN(1.)
    197188c
    198       PRINT 990
     189      WRITE(6,990)
    199190
    200191c     ----------------------------------------------------------------
    … …  
    205196       IF( ysinus )  THEN
    206197c
    207         PRINT *,' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '
     198        WRITE(6,*) ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '
    208199c
    209200c   .... utilisation de f(x,y )  avec  y  =  sinus de la latitude  .....
    … …  
    215206       ELSE
    216207c
    217         PRINT *,' *** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
     208        WRITE(6,*) '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
    218209
    219210c  .... utilisation  de f(x,y) a tangente sinusoidale , y etant la latit. ...
    … …  
    267258      ELSE
    268259c
    269 c   ....  utilisation  de fxyhyper , f(x,y) a derivee tangente hyperbol.
     260c   ....  Utilisation  de fxyhyper , f(x,y) a derivee tangente hyperbol.
    270261c   .....................................................................
    271262
    272        PRINT *,'*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
    273 
    274        CALL fxyhyper( clat, grossismy, dzoomy, tauy, deltay  ,
    275      ,               clon, grossismx, dzoomx, taux          ,
     263      WRITE(6,*)'*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
     264 
     265       CALL fxyhyper( clat, grossismy, dzoomy, tauy    ,
     266     ,                clon, grossismx, dzoomx, taux    ,
    276267     , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2  ,
    277268     , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
    … …  
    387378c
    388379      IF ( j. eq. jjp1 )  THEN
    389 
    390380       yprp               = yprimu2(j-1)
    391381       rlatp              = rlatu2 (j-1)
    392 cc       yprp             = fyprim( FLOAT(j) - 0.25 )
    393 cc       rlatp            = fy    ( FLOAT(j) - 0.25 )
     382ccc       yprp             = fyprim( FLOAT(j) - 0.25 )
     383ccc       rlatp            = fy    ( FLOAT(j) - 0.25 )
    394384c
    395385      coslatp             = COS( rlatp )
    … …  
    425415        rlatm    = rlatu1 (  j  )
    426416        yprm     = yprimu1(  j  )
    427 c         rlatp    = fy    ( FLOAT(j) - 0.25 )
    428 c         yprp     = fyprim( FLOAT(j) - 0.25 )
    429 c         rlatm    = fy    ( FLOAT(j) + 0.25 )
    430 c         yprm     = fyprim( FLOAT(j) + 0.25 )
     417cc         rlatp    = fy    ( FLOAT(j) - 0.25 )
     418cc         yprp     = fyprim( FLOAT(j) - 0.25 )
     419cc         rlatm    = fy    ( FLOAT(j) + 0.25 )
     420cc         yprm     = fyprim( FLOAT(j) + 0.25 )
    431421
    432422         coslatm  = COS( rlatm )
    … …  
    490480  36  CONTINUE
    491481c
    492 c   ....  Modif  P. Le Van  ( 4/07/96 )  .....
    493482c
    494483      aire    (iip1,j) = aire    (1,j)
    … …  
    620609         aiuscv2gam(iip1,j)  = aiuscv2gam(1,j)
    621610      ENDDO
    622 c
     611
    623612c
    624613c   calcul des aires aux  poles :
    … …  
    666655c-----------------------------------------------------------------------
    667656c
    668        PRINT *,' INIGEOM  RLONV '
     657       WRITE(6,*) '   ***  Coordonnees de la grille  *** '
     658       WRITE(6,995)
     659c
     660       WRITE(6,*) '   LONGITUDES  aux pts.   V  ( degres )  '
     661       WRITE(6,995)
    669662        DO i=1,iip1
    670663         rlonvv(i) = rlonv(i)*180./pi
    671664        ENDDO
    672        PRINT 400,rlonvv
    673 c
    674        PRINT *,' RLATV '
     665       WRITE(6,400) rlonvv
     666c
     667       WRITE(6,995)
     668       WRITE(6,*) '   LATITUDES   aux pts.   V  ( degres )  '
     669       WRITE(6,995)
    675670        DO i=1,jjm
    676671         rlatuu(i)=rlatv(i)*180./pi
    677672        ENDDO
    678        PRINT 400,(rlatuu(i),i=1,jjm)
     673       WRITE(6,400) (rlatuu(i),i=1,jjm)
    679674c
    680675        DO i=1,iip1
    681676          rlonvv(i)=rlonu(i)*180./pi
    682677        ENDDO
    683        PRINT *,' RLONU '
    684        PRINT 400,rlonvv
    685 c
    686        PRINT *,' RLATU '
     678       WRITE(6,995)
     679       WRITE(6,*) '   LONGITUDES  aux pts.   U  ( degres )  '
     680       WRITE(6,995)
     681       WRITE(6,400) rlonvv
     682       WRITE(6,995)
     683
     684       WRITE(6,*) '   LATITUDES   aux pts.   U  ( degres )  '
     685       WRITE(6,995)
    687686        DO i=1,jjp1
    688687         rlatuu(i)=rlatu(i)*180./pi
    689688        ENDDO
    690        PRINT 400,(rlatuu(i),i=1,jjp1)
    691 c
     689       WRITE(6,400) (rlatuu(i),i=1,jjp1)
     690       WRITE(6,995)
     691c
     692444    format(f10.3,f6.0)
    692693400    FORMAT(1x,8f8.2)
    693694990    FORMAT(//)
     695995    FORMAT(/)
    694696c
    695697      RETURN

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/paramet.h

    r2 r232  
    77      INTEGER jcfil,jcfllm
    88
    9       PARAMETER( iip1= iim+1, iip2= iim+2, iip3= iim+3, jjp1= jjm+1 )
     9      PARAMETER( iip1= iim+1-1/iim,iip2=iim+2,iip3=iim+3
     10     s    ,jjp1=jjm+1-1/jjm)
    1011      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
    1112      PARAMETER( kftd  = iim/2 -ndm )

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/physdem.F

    r115 r232  
    527527      ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
    528528#endif
    529 c
    530529      ierr = NF_REDEF (nid)
    531530#ifdef NC_DOUBLE
    … …  
    540539      ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
    541540#endif
    542 c
    543541      ierr = NF_REDEF (nid)
    544542#ifdef NC_DOUBLE
    … …  
    553551      ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
    554552#endif
    555 c
    556553      ierr = NF_REDEF (nid)
    557554#ifdef NC_DOUBLE
    … …  
    566563      ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
    567564#endif
    568 c
    569565      ierr = NF_REDEF (nid)
    570566#ifdef NC_DOUBLE
    … …  
    579575      ierr = NF_PUT_VAR_REAL (nid,nvarid,zpic)
    580576#endif
    581 c
    582577      ierr = NF_REDEF (nid)
    583578#ifdef NC_DOUBLE
    … …  
    592587      ierr = NF_PUT_VAR_REAL (nid,nvarid,zval)
    593588#endif
    594 c
    595589      ierr = NF_REDEF (nid)
    596590#ifdef NC_DOUBLE

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/serre.h

    r2 r232  
     1c
     2c $Header$
     3c
    14c..include serre.h
    25c
    36       REAL clon,clat,transx,transy,alphax,alphay,pxo,pyo,
    4      ,  grossismx, grossismy, dzoomx, dzoomy
     7     ,  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
    58       COMMON/serre/clon,clat,transx,transy,alphax,alphay,pxo,pyo ,
    6      ,  grossismx, grossismy, dzoomx, dzoomy
     9     ,  grossismx, grossismy, dzoomx, dzoomy,taux,tauy

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/tracvl.F

    r2 r232  
    11      SUBROUTINE tracvl(numvanle,iapp_tracvl,nq,pbaru,pbarv ,
    2      *                            p,  masse , q, iapptrac    )
     2     *                    p, masse , q, iapptrac, iadv1, teta, pk  )
    33c
    44c     Auteur :  F. Hourdin
    … …  
    66c
    77ccc   ..   Modif. P. Le Van  ( 20/12/97 )  ...
     8c                 F. Codron     (10/99)
     9
    810c
    911      IMPLICIT NONE
    … …  
    1517#include "comgeom.h"
    1618
    17       INTEGER nq,iapp_tracvl
     19c     .... Arguments  ....
     20c
     21      INTEGER numvanle, nq, iapp_tracvl, iapptrac, iadv1
    1822
    1923      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
    2024      REAL q(ip1jmp1,llm,nq),masse(ip1jmp1,llm)
    21       REAL p( ip1jmp1,llmp1 )
     25      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
     26      REAL pk(ip1jmp1,llm)
    2227
     28c     ....  var. locales  .....
     29c
    2330      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
    2431      REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
    … …  
    2633      REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
    2734
    28       INTEGER iapptrac
    2935
    30       INTEGER iadvtr, numvanle
     36      INTEGER iadvtr, numvan
    3137      INTEGER ij,l,iq
    3238      REAL zdpmin, zdpmax
    3339      EXTERNAL  minmax
    34  
    35       SAVE iadvtr, massem,pbaruc,pbarvc
     40      SAVE iadvtr, massem, pbaruc, pbarvc, numvan
    3641      DATA iadvtr/0/
     42
     43      numvan = numvanle
    3744
    3845      IF(iadvtr.EQ.0) THEN
    … …  
    99106
    100107c   Advection proprement dite.
    101          DO iq = numvanle, nq
     108c
     109c   test sur iadv1 pour le schema de vapeur d'eau
     110c
     111         IF (numvanle.EQ.1.AND.iadv1.EQ.4) THEN
     112           CALL vlspltqs( q(1,1,1), 2., massem, wg ,
     113     *                 pbarug,pbarvg,dtvr,p,pk,teta )
     114           numvan = 2
     115         ENDIF
     116
     117         DO iq = numvan, nq
    102118          CALL vlsplt( q(1,1,iq), 2. ,massem,wg,pbarug,pbarvg,dtvr )
    103119         ENDDO

  • LMDZ.3.3/branches/rel-LF/libf/dyn3d/vanleer.F

    r2 r232  
    11      SUBROUTINE vanleer(numvanle,iapp_tracvl,nq,q,pbaru,pbarv ,
    2      *                                    p ,masse,  dq          )
     2     *                     p ,masse, dq ,  iadv1, teta, pk      )
    33c
    44      IMPLICIT NONE
    55c
    6 c     Auteurs:   F.Hourdin , P.Le Van, F.Forget
     6c     Auteurs:   F.Hourdin , P.Le Van, F.Forget, F.Codron 
     7c
     8c     F.Codron (10/99) : ajout humidite specifique pour eau vapeur
    79c=======================================================================
    810c
    … …  
    1820c   Arguments:
    1921c   ----------
    20       INTEGER nq,numvanle,iapp_tracvl
     22      INTEGER nq, numvanle, iapp_tracvl, iadv1
    2123      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm),masse(ip1jmp1,llm)
    2224      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nq),dq( ip1jmp1,llm,nq )
    23 
     25      REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
    2426c  ..................................................................
    2527c
    … …  
    2830c     numvanle est le numero du premier traceur qui appelle le
    2931c       shema de Van-Leer  (  1<=  numvanle <=  nqmx )  .
     32c     si numvanle=1, iadv1=3 si Van-Leer seul, iadv1=4 si
     33c     Van-Leer+humidite specifique.
    3034c  ..................................................................
    3135c
    … …  
    6165
    6266      CALL tracvl( numvanle,iapp_tracvl,nq,pbaru,pbarv,p , masse,q  ,
    63      *                                iapptrac                       )
     67     *                      iapptrac, iadv1, teta ,pk              )
    6468
    6569      IF( numvanle.EQ.1 ) THEN
Note: See TracChangeset for help on using the changeset viewer.