|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 C Jan-2011 P. Marguinaud Thread-safe FA 00002 SUBROUTINE FACDEC_MT (FA, KREP, PA, PMIN, KNBIT, KDEC) 00003 USE FA_MOD, ONLY : FA_COM 00004 USE PARKIND1, ONLY : JPRB 00005 USE YOMHOOK , ONLY : LHOOK, DR_HOOK 00006 !**** 00007 ! Sous-programme de calcul du FACteur d'echelle DECimal optimal 00008 ! pour coder un champ d'amplitude PA donnee, sur KNBIT bits avec 00009 ! une precision maximale, compte tenu d'un intervalle fixe 00010 ! de facteurs decimaux. 00011 !** 00012 ! Arguments : KREP (Sortie) ==> Code-reponse du sous-programme; 00013 ! 0: OK; -1: pb rencontre 00014 ! PA (Entree) ==> Amplitude du champ a compacter; 00015 ! PMIN (Entree) ==> Minimum du champ a compacter; 00016 ! KNBIT (Entree) ==> Nb de bits servant au compactage; 00017 ! KDEC (Sortie) ==> Facteur d'echelle decimal optimal; 00018 ! 00019 ! Modifications: 00020 ! R. El Khatib 12-Aug-2009 Bugfix for compatibility with Gribex 00021 ! 00022 IMPLICIT NONE 00023 #include "precision.h" 00024 TYPE(FA_COM) :: FA 00025 REAL (KIND=JPDBLR) PA, PMIN 00026 INTEGER KREP, KNBIT, KDEC 00027 ! 00028 REAL (KIND=JPDBLR) XNBINT, XTINYR4, XHUGER4 00029 INTEGER IDECMIN, IDECMAX, INBINT 00030 INTEGER JDEC, IE, INUTIL, INUMAX, IEMAX 00031 INTEGER INU0, IE0 00032 #include "facom_mt.h" 00033 ! 00034 !** 00035 ! 1. - CONTROLES ET INITIALISATIONS. 00036 !----------------------------------------------------------------------- 00037 ! 00038 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00039 IF (LHOOK) CALL DR_HOOK('FACDEC_MT',0,ZHOOK_HANDLE) 00040 IF (KNBIT.LE.0 .OR. KNBIT.GT.64) THEN 00041 KREP=-1 00042 WRITE (UNIT=FA%NULOUT,FMT=*)'****' 00043 WRITE (UNIT=FA%NULOUT,FMT=*) 00044 S '**** FACDEC: ERROR, bits number out of range 1-64' 00045 WRITE (UNIT=FA%NULOUT,FMT=*)'**** KNBIT = ',KNBIT 00046 WRITE (UNIT=FA%NULOUT,FMT=*) 00047 S '**** ! Optimal decimal scale factor is not computed !' 00048 WRITE (UNIT=FA%NULOUT,FMT=*)'****' 00049 IF (LHOOK) CALL DR_HOOK('FACDEC_MT',1,ZHOOK_HANDLE) 00050 RETURN 00051 ENDIF 00052 IF (ABS(PA) .LT. TINY(PA)) THEN 00053 KREP = 0 00054 KDEC = 0 00055 IF (FA%LFAMOP) THEN 00056 WRITE (UNIT=FA%NULOUT,FMT=*)'////' 00057 WRITE (UNIT=FA%NULOUT,FMT=*) 00058 S '//// FACDEC: WARNING, range of the field is null :', 00059 S PA 00060 WRITE (UNIT=FA%NULOUT,FMT=*)'////' 00061 ENDIF 00062 IF (LHOOK) CALL DR_HOOK('FACDEC_MT',1,ZHOOK_HANDLE) 00063 RETURN 00064 ENDIF 00065 KREP = 0 00066 INUMAX = 0 00067 IDECMIN = -15 00068 IDECMAX = 5 00069 XTINYR4 = TINY(1._4) 00070 XHUGER4 = HUGE(1._4) 00071 INBINT = 2**KNBIT -1 00072 XNBINT = FLOAT(INBINT) 00073 ! Cas du facteur decimal nul (reference a calculer dans tous les cas) 00074 JDEC = 0 00075 CALL FACTEC_MT(FA, KREP, PA, KNBIT, JDEC, IE0, INU0) 00076 00077 ! 00078 !** 00079 ! 2. - BOUCLE TESTANT LES FACTEURS DECIMAUX 00080 !----------------------------------------------------------------------- 00081 ! 00082 !!OCL SCALAR 00083 DO JDEC = IDECMIN, IDECMAX 00084 ! 00085 ! 3 tests sur la pertinence de ce facteur decimal 00086 ! 00087 ! 0/ PA * 10**JDEC > 1.E-12 00088 ! Du fait d'un test exagere dans gribex 00089 IF (PA * 10._JPDBLR**JDEC .LE. 1.E-12) CYCLE 00090 ! 1/ PMIN * 10**JDEC > TINY(real*4) pour decodage eventuel 00091 ! de la valeur de reference en arithmetique 32 bits IEEE. 00092 IF (ABS(PMIN) .GT. TINY(PMIN)) THEN 00093 IF (LOG10(ABS(PMIN))+JDEC .LE. LOG10( XTINYR4 ) ) CYCLE 00094 ELSE 00095 ! instruction "bidon" pour que l'optimisation ne fasse pas 00096 ! l'evaluation de LOG10(ABS(PMIN)) par anticipation (Pb sur VPP5000) 00097 PMIN = 0._JPDBLR 00098 ENDIF 00099 ! 2/ PA * 10**JDEC inclus dans le domaine de validite des real*8 00100 IF (ABS(LOG10(ABS(PA))+JDEC) .GE. RANGE(PA)) CYCLE 00101 ! 00102 CALL FACTEC_MT(FA, KREP, PA, KNBIT, JDEC, IE, INUTIL) 00103 IF (KREP.NE.0) CYCLE 00104 ! 3/ PMIN*10**JDEC + (2**KNBIT-1)*2**IE < HUGE(real*4) pour decodage 00105 ! eventuel du max du champ en arithmetique 32 bits IEEE. 00106 IF (PMIN*10._JPDBLR**JDEC + XNBINT*2._JPDBLR**IE .GE. 00107 S XHUGER4) CYCLE 00108 ! 4/ Le facteur d'echelle binaire doit pouvoir etre code sur 1 octet, 00109 ! il doit donc etre compris entre -126 et 127: 00110 IF (IE.LT.-126 .OR. IE.GT.127) CYCLE 00111 ! 00112 IF (INUTIL.GT.INUMAX) THEN 00113 INUMAX = INUTIL 00114 IEMAX = IE 00115 KDEC = JDEC 00116 ENDIF 00117 ENDDO 00118 ! 00119 IF (INUMAX.EQ.0) THEN 00120 KREP=-1 00121 WRITE (UNIT=FA%NULOUT,FMT=*)'****' 00122 WRITE (UNIT=FA%NULOUT,FMT=*) 00123 S '**** FACDEC: all the decimal factors comprised between' 00124 WRITE (UNIT=FA%NULOUT,FMT=*) 00125 S '**** ',IDECMIN,' and ',IDECMAX,' are rejected !!' 00126 WRITE (UNIT=FA%NULOUT,FMT=*) 00127 S '**** Range and min of the field are :', PA, PMIN 00128 WRITE (UNIT=FA%NULOUT,FMT=*) 00129 S '****' 00130 WRITE (UNIT=FA%NULOUT,FMT=*)'****' 00131 IF (LHOOK) CALL DR_HOOK('FACDEC_MT',1,ZHOOK_HANDLE) 00132 RETURN 00133 ENDIF 00134 IF (FA%LFAMOP) THEN 00135 WRITE (UNIT=FA%NULOUT,FMT=*) 00136 S 'FACDEC: champ d''amplitude ',PA,' ,de minimum ',PMIN 00137 WRITE (UNIT=FA%NULOUT,FMT=*) 00138 S ' => fact decimal opt de',KDEC, 00139 S ' ,pour 1 fact binaire de ',IEMAX 00140 WRITE (UNIT=FA%NULOUT,FMT='(1X,A,I3,A,I9,A,I9,A,F5.1,A,E11.4)') 00141 S ' Eff des', 00142 S KNBIT,' bits = ',INUMAX,' sur ',INBINT,' soit: ', 00143 S FLOAT(INUMAX)*100./XNBINT,' % et une precision de ', 00144 S 10.**(-KDEC)*2.**(IEMAX-1) 00145 WRITE (UNIT=FA%NULOUT,FMT=*) 00146 S ' a comparer, si le fact decimal = 0, avec' 00147 WRITE (UNIT=FA%NULOUT,FMT='(1X,A,I9,A,I9,A,F5.1,A,E11.4)') 00148 S ' une efficacite de ', 00149 S INU0,' sur ',INBINT,' soit: ',FLOAT(INU0)*100./XNBINT, 00150 S ' % et une precision de ',2.**(IE0-1) 00151 ENDIF 00152 ! 00153 IF (LHOOK) CALL DR_HOOK('FACDEC_MT',1,ZHOOK_HANDLE) 00154 END 00155
1.8.0
