Changeset 270
- Timestamp:
- Aug 18, 2011, 12:09:27 PM (13 years ago)
- Location:
- trunk
- Files:
-
- 2 added
- 42 edited
- 4 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/DOC/chantiers/commit_importants.log
r234 r270 949 949 ********************* 950 950 plus de meso_ files dans la physique martienne. physiq et inifis sont communs au GCM avec des precompiler flags. 951 952 ********************* 953 **** commit_v270 **** 954 ********************* 955 Ehouarn: mises à jour des dynamiques (seq et //) pour suivre les 956 développements de LMDZ5 terrestre (LMDZ5 trunk, rev 1560) 957 958 * scripts & co.: changement (mineur) dans "makelmdz_fcm" et "makelmdz" et les 959 fichiers "arch*path" sur la façon dont sont définies et 960 utilisées les références à NetCDF 961 962 * dyn3d: 963 - advtrac.F devient advtrac.F90 964 - inidissip.F devient inidissip.F90 + utilisation de "dissip_period" au lieu de 965 "idissip" et calcul automatique de "dissip_period" si dissip_period=0 966 dans gcm.def 967 - iniacademic.F90 : adaptation pour aquaplanètes terrestres (iflag_phys>2) 968 - ajout de grilles_gcm_netcdf_sub.F90 (pas utile pour les planètes) 969 - leapfrog.F :ajout du logical 'physic' (pour régler l'appel à la physique, 970 avec les cas particuliers aquaplanètes) + correction bug sur cumul des 971 flux de masse pour l'advection des tracers. 972 - ini_paramLMDZ_dyn.h : remplacement de "idissip" par "dissip_period" 973 - comdissipn.h : passage a un format compatible fixe/libre 974 - control_mod.F90 : remplacement de "idissip" par "dissip_period" 975 - defrun.F : remplacement de "idissip" par "dissip_period" 976 - gcm.F : adaptations pour aquaplanètes terrestres (iflag_phys>2) 977 - iniconst.F : dtdiss est calculé dans inidissip.F90 978 - vlsplt.F : changements cosmetiques 979 - ce0l.F90 : ajout appel à grilles_gcm_netcdf_sub 980 - integrd.F : changements cosmetiques 981 - limit_netcdf.F90 : mise à jour de la version terrestre (pas important pour 982 les planètes) 983 - conf_gcm.F : ajout lecture du flag 'dissip_period' (au lieu de idissip) et 984 du flag 'grilles_gcm_netcdf' (pas important pour les planètes) 985 - write_paramLMDZ_dyn.h : remplacement de "idissip" par "dissip_period" 986 987 *dyn3dpar: 988 - Idem que pour dyn3d pour le fichiers identiques entre dyn3d et dyn3dpar: 989 comdissipn.h, control_mod.F90, iniconst.F, grilles_gcm_netcdf_sub.F90, 990 ce0l.F90, inidissip.F90, iniacademic.F90, limit_netcdf.F90, conf_gcm.F, 991 defrun.F 992 - advtrac_p.F devient advtrac_P.F90 993 - integrd_p.F : changements cosmetiques 994 - leapfrog_p.F90 : changements cosmetiques + correction bug sur cumul des 995 flux de masse 996 - bilan_dyn_p.F : correction bug mineur 997 - exner_hyb_p.F : correction bug OMP majeur 998 - exner_milieu_p.F : correction bug OMP majeur -
trunk/DOC/chantiers/meschantiers-Ehouarn.txt
r127 r270 12 12 puis mes modifs en r124 et finalement r126) donc a priori 13 13 OK maintenant 14 15 - dernière mise à jour des dynamiques (par rapport au trunk de LMDZ5 16 "Terrestre" rev. 1560) avec la rev 270. 17 Reste à faire l'exercice de validation des modifs en //... 14 18 -
trunk/LMDZ.COMMON/arch/arch-AMD64_CICLAD.path
r52 r270 1 NETCDF_LIBDIR= /opt/netcdf/pgf/lib2 NETCDF_INCDIR= /opt/netcdf/pgf/include1 NETCDF_LIBDIR="-L/opt/netcdf/pgf/lib -lnetcdf" 2 NETCDF_INCDIR=-I/opt/netcdf/pgf/include 3 3 IOIPSL_INCDIR=$trunk/ioipsl/lib 4 4 IOIPSL_LIBDIR=$trunk/ioipsl/lib -
trunk/LMDZ.COMMON/arch/arch-ES_MOON.path
r1 r270 1 NETCDF_LIBDIR= /S/home010/c0010/ES/lib2 NETCDF_INCDIR= /S/home010/c0010/ES/include1 NETCDF_LIBDIR="-L/S/home010/c0010/ES/lib -lnetcdf" 2 NETCDF_INCDIR=-I/S/home010/c0010/ES/include 3 3 IOIPSL_INCDIR=$LMDGCM/../../lib 4 4 IOIPSL_LIBDIR=$LMDGCM/../../lib -
trunk/LMDZ.COMMON/arch/arch-GNOMEs.path
r130 r270 1 NETCDF_LIBDIR= /usr/local/lib2 NETCDF_INCDIR= /usr/local/include1 NETCDF_LIBDIR="-L/usr/local/lib -lnetcdf" 2 NETCDF_INCDIR=-I/usr/local/include 3 3 IOIPSL_INCDIR=$HOME/LMDZ5/ioipsl/lib 4 4 IOIPSL_LIBDIR=$HOME/LMDZ5/ioipsl/lib -
trunk/LMDZ.COMMON/arch/arch-IA64_PLATINE.path
r1 r270 1 NETCDF_LIBDIR= '/usr/lib -lnetcdff'2 NETCDF_INCDIR= /usr/include1 NETCDF_LIBDIR="-L/usr/lib -lnetcdff -lnetcdf" 2 NETCDF_INCDIR=-I/usr/include 3 3 IOIPSL_INCDIR=$LMDGCM/../../lib 4 4 IOIPSL_LIBDIR=$LMDGCM/../../lib -
trunk/LMDZ.COMMON/arch/arch-PW6_VARGAS.path
r1 r270 1 NETCDF_LIBDIR= /usr/local/pub/NetCDF/3.6.3/lib2 NETCDF_INCDIR= /usr/local/pub/NetCDF/3.6.3/include1 NETCDF_LIBDIR="${NETCDF_LDFLAGS:--L/usr/local/pub/NetCDF/3.6.2/lib -lnetcdf}" 2 NETCDF_INCDIR="${NETCDF_FFLAGS:--I/usr/local/pub/NetCDF/3.6.2/include}" 3 3 IOIPSL_INCDIR=$LMDGCM/../../lib 4 4 IOIPSL_LIBDIR=$LMDGCM/../../lib -
trunk/LMDZ.COMMON/arch/arch-SX8_BRODIE.path
r1 r270 1 NETCDF_LIBDIR= /SXlocal/pub/netCDF/3.6.1-openmp/lib2 NETCDF_INCDIR= /SXlocal/pub/netCDF/3.6.1-openmp/include1 NETCDF_LIBDIR="-L/SXlocal/pub/netCDF/3.6.1-openmp/lib -lnetcdf" 2 NETCDF_INCDIR=-I/SXlocal/pub/netCDF/3.6.1-openmp/include 3 3 IOIPSL_INCDIR=$LMDGCM/../../lib 4 4 IOIPSL_LIBDIR=$LMDGCM/../../lib -
trunk/LMDZ.COMMON/arch/arch-SX8_MERCURE.path
r1 r270 1 NETCDF_LIBDIR= ${NETCDF_SX_LIBDIR:-/usr/local/SX8/soft/netcdf/lib}2 NETCDF_INCDIR= ${NETCDF_SX_INCLUDEDIR:-/usr/local/SX8/soft/netcdf/include}1 NETCDF_LIBDIR="-L ${NETCDF_SX_LIBDIR:-/usr/local/SX8/soft/netcdf/lib} -lnetcdf" 2 NETCDF_INCDIR=-I${NETCDF_SX_INCLUDEDIR:-/usr/local/SX8/soft/netcdf/include} 3 3 IOIPSL_INCDIR=$LMDGCM/../../lib 4 4 IOIPSL_LIBDIR=$LMDGCM/../../lib -
trunk/LMDZ.COMMON/arch/arch-SX9_MERCURE.path
r1 r270 1 NETCDF_LIBDIR= ${NETCDF_SX_LIBDIR:-/ccc/applications/sx9/netcdf-3.6.1/lib}2 NETCDF_INCDIR= ${NETCDF_SX_INCLUDEDIR:-/ccc/applications/sx9/netcdf-3.6.1/include}1 NETCDF_LIBDIR="-L${NETCDF_SX_LIBDIR:-/ccc/applications/sx9/netcdf-3.6.1/lib} -lnetcdf" 2 NETCDF_INCDIR=-I${NETCDF_SX_INCLUDEDIR:-/ccc/applications/sx9/netcdf-3.6.1/include} 3 3 IOIPSL_INCDIR=$LMDGCM/../../lib 4 4 IOIPSL_LIBDIR=$LMDGCM/../../lib -
trunk/LMDZ.COMMON/arch/arch-X64_TITANE.path
r1 r270 1 NETCDF_LIBDIR=" $NETCDF_LIB_DIR -lnetcdff"2 NETCDF_INCDIR= $NETCDF_INC_DIR1 NETCDF_LIBDIR="-L$NETCDF_LIB_DIR -lnetcdff -lnetcdf" 2 NETCDF_INCDIR=-I$NETCDF_INC_DIR 3 3 IOIPSL_INCDIR=$LMDGCM/../../lib 4 4 IOIPSL_LIBDIR=$LMDGCM/../../lib -
trunk/LMDZ.COMMON/arch/arch-linux-32bit.path
r1 r270 1 NETCDF_LIBDIR= /usr/local/netcdf-pgi/lib2 NETCDF_INCDIR= /usr/local/netcdf-pgi/include1 NETCDF_LIBDIR="-L/usr/local/netcdf-pgi/lib -lnetcdf" 2 NETCDF_INCDIR=-I/usr/local/netcdf-pgi/include 3 3 IOIPSL_INCDIR=/usr/local/guez/modipsl/lib 4 4 IOIPSL_LIBDIR=/usr/local/guez/modipsl/lib -
trunk/LMDZ.COMMON/arch/arch-linux-64bit.path
r130 r270 1 NETCDF_LIBDIR= /donnees/emlmd/netcdf64-4.0.1_pgi/lib2 NETCDF_INCDIR= /donnees/emlmd/netcdf64-4.0.1_pgi/include1 NETCDF_LIBDIR="-L/donnees/emlmd/netcdf64-4.0.1_pgi/lib -lnetcdf" 2 NETCDF_INCDIR=-I/donnees/emlmd/netcdf64-4.0.1_pgi/include 3 3 IOIPSL_INCDIR=/donnees/emlmd/modipsl64_pgi/modipsl/lib 4 4 IOIPSL_LIBDIR=/donnees/emlmd/modipsl64_pgi/modipsl/lib -
trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90
r269 r270 1 ! 2 ! $Id: advtrac.F 1446 2010-10-22 09:27:25Z emillour $ 3 ! 4 c 5 c 6 SUBROUTINE advtrac(pbaru,pbarv , 7 * p, masse,q,iapptrac,teta, 8 * flxw, 9 * pk) 10 c Auteur : F. Hourdin 11 c 12 c Modif. P. Le Van (20/12/97) 13 c F. Codron (10/99) 14 c D. Le Croller (07/2001) 15 c M.A Filiberti (04/2002) 16 c 17 USE infotrac 18 USE control_mod 19 20 21 IMPLICIT NONE 22 c 23 #include "dimensions.h" 24 #include "paramet.h" 25 #include "comconst.h" 26 #include "comvert.h" 27 #include "comdissip.h" 28 #include "comgeom2.h" 29 #include "logic.h" 30 #include "temps.h" 31 #include "ener.h" 32 #include "description.h" 33 #include "iniprint.h" 34 35 c------------------------------------------------------------------- 36 c Arguments 37 c------------------------------------------------------------------- 38 c Ajout PPM 39 c-------------------------------------------------------- 40 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) 41 c-------------------------------------------------------- 42 INTEGER iapptrac 43 REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 44 REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm) 45 REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm) 46 REAL pk(ip1jmp1,llm) 47 REAL flxw(ip1jmp1,llm) 48 49 c------------------------------------------------------------- 50 c Variables locales 51 c------------------------------------------------------------- 52 53 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) 54 REAL massem(ip1jmp1,llm),zdp(ip1jmp1) 55 REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 56 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu 57 INTEGER iadvtr 58 INTEGER ij,l,iq,iiq 59 REAL zdpmin, zdpmax 60 EXTERNAL minmax 61 SAVE iadvtr, massem, pbaruc, pbarvc 62 DATA iadvtr/0/ 63 c---------------------------------------------------------- 64 c Rajouts pour PPM 65 c---------------------------------------------------------- 66 INTEGER indice,n 67 REAL dtbon ! Pas de temps adaptatif pour que CFL<1 68 REAL CFLmaxz,aaa,bbb ! CFL maximum 69 REAL psppm(iim,jjp1) ! pression au sol 70 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) 71 REAL qppm(iim*jjp1,llm,nqtot) 72 REAL fluxwppm(iim,jjp1,llm) 73 REAL apppm(llmp1), bpppm(llmp1) 74 LOGICAL dum,fill 75 DATA fill/.true./ 76 DATA dum/.true./ 77 78 integer,save :: countcfl=0 79 real cflx(ip1jmp1,llm) 80 real cfly(ip1jm,llm) 81 real cflz(ip1jmp1,llm) 82 real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm) 83 84 IF(iadvtr.EQ.0) THEN 85 CALL initial0(ijp1llm,pbaruc) 86 CALL initial0(ijmllm,pbarvc) 87 ENDIF 88 89 c accumulation des flux de masse horizontaux 90 DO l=1,llm 91 DO ij = 1,ip1jmp1 92 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 93 ENDDO 94 DO ij = 1,ip1jm 95 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 96 ENDDO 97 ENDDO 98 99 c selection de la masse instantannee des mailles avant le transport. 100 IF(iadvtr.EQ.0) THEN 101 102 CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) 103 ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 104 c 105 ENDIF 106 107 iadvtr = iadvtr+1 108 iapptrac = iadvtr 109 110 111 c Test pour savoir si on advecte a ce pas de temps 112 IF ( iadvtr.EQ.iapp_tracvl ) THEN 113 114 cc .. Modif P.Le Van ( 20/12/97 ) .... 115 cc 116 117 c traitement des flux de masse avant advection. 118 c 1. calcul de w 119 c 2. groupement des mailles pres du pole. 120 121 CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) 122 123 ! ... Flux de masse diaganostiques traceurs 124 flxw = wg / REAL(iapp_tracvl) 125 126 c test sur l'eventuelle creation de valeurs negatives de la masse 127 DO l=1,llm-1 128 DO ij = iip2+1,ip1jm 129 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) 130 s - pbarvg(ij-iip1,l) + pbarvg(ij,l) 131 s + wg(ij,l+1) - wg(ij,l) 132 ENDDO 133 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 134 DO ij = iip2,ip1jm 135 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 136 ENDDO 137 138 139 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 140 141 IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN 142 PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, 143 s ' MAX:', zdpmax 144 ENDIF 145 146 ENDDO 147 148 149 c------------------------------------------------------------------- 150 ! Calcul des criteres CFL en X, Y et Z 151 c------------------------------------------------------------------- 152 153 if (countcfl == 0. ) then 154 cflxmax(:)=0. 155 cflymax(:)=0. 156 cflzmax(:)=0. 157 endif 158 159 countcfl=countcfl+iapp_tracvl 160 cflx(:,:)=0. 161 cfly(:,:)=0. 162 cflz(:,:)=0. 163 do l=1,llm 164 do ij=iip2,ip1jm-1 165 if (pbarug(ij,l)>=0.) then 166 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) 167 else 168 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) 169 endif 170 enddo 171 enddo 172 do l=1,llm 173 do ij=iip2,ip1jm-1,iip1 174 cflx(ij+iip1,l)=cflx(ij,l) 175 enddo 176 enddo 177 178 do l=1,llm 179 do ij=1,ip1jm 180 if (pbarvg(ij,l)>=0.) then 181 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) 182 else 183 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) 184 endif 185 enddo 186 enddo 187 188 do l=2,llm 189 do ij=1,ip1jm 190 if (wg(ij,l)>=0.) then 191 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) 192 else 193 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) 194 endif 195 enddo 196 enddo 197 198 do l=1,llm 199 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) 200 cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) 201 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) 202 enddo 1 ! $Id: advtrac.F90 1549 2011-07-05 08:41:12Z lguez $ 2 3 SUBROUTINE advtrac(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk) 4 ! Auteur : F. Hourdin 5 ! 6 ! Modif. P. Le Van (20/12/97) 7 ! F. Codron (10/99) 8 ! D. Le Croller (07/2001) 9 ! M.A Filiberti (04/2002) 10 ! 11 USE infotrac 12 USE control_mod 13 14 15 IMPLICIT NONE 16 ! 17 include "dimensions.h" 18 include "paramet.h" 19 include "comconst.h" 20 include "comvert.h" 21 include "comdissip.h" 22 include "comgeom2.h" 23 include "logic.h" 24 include "temps.h" 25 include "ener.h" 26 include "description.h" 27 include "iniprint.h" 28 29 !------------------------------------------------------------------- 30 ! Arguments 31 !------------------------------------------------------------------- 32 ! Ajout PPM 33 !-------------------------------------------------------- 34 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) 35 !-------------------------------------------------------- 36 INTEGER iapptrac 37 REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 38 REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm) 39 REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm) 40 REAL pk(ip1jmp1,llm) 41 REAL flxw(ip1jmp1,llm) 42 43 !------------------------------------------------------------- 44 ! Variables locales 45 !------------------------------------------------------------- 46 47 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) 48 REAL massem(ip1jmp1,llm),zdp(ip1jmp1) 49 REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 50 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu 51 INTEGER iadvtr 52 INTEGER ij,l,iq,iiq 53 REAL zdpmin, zdpmax 54 EXTERNAL minmax 55 SAVE iadvtr, massem, pbaruc, pbarvc 56 DATA iadvtr/0/ 57 !---------------------------------------------------------- 58 ! Rajouts pour PPM 59 !---------------------------------------------------------- 60 INTEGER indice,n 61 REAL dtbon ! Pas de temps adaptatif pour que CFL<1 62 REAL CFLmaxz,aaa,bbb ! CFL maximum 63 REAL psppm(iim,jjp1) ! pression au sol 64 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) 65 REAL qppm(iim*jjp1,llm,nqtot) 66 REAL fluxwppm(iim,jjp1,llm) 67 REAL apppm(llmp1), bpppm(llmp1) 68 LOGICAL dum,fill 69 DATA fill/.true./ 70 DATA dum/.true./ 71 72 integer,save :: countcfl=0 73 real cflx(ip1jmp1,llm) 74 real cfly(ip1jm,llm) 75 real cflz(ip1jmp1,llm) 76 real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm) 77 78 IF(iadvtr.EQ.0) THEN 79 CALL initial0(ijp1llm,pbaruc) 80 CALL initial0(ijmllm,pbarvc) 81 ENDIF 82 83 ! accumulation des flux de masse horizontaux 84 DO l=1,llm 85 DO ij = 1,ip1jmp1 86 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 87 ENDDO 88 DO ij = 1,ip1jm 89 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 90 ENDDO 91 ENDDO 92 93 ! selection de la masse instantannee des mailles avant le transport. 94 IF(iadvtr.EQ.0) THEN 95 96 CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) 97 !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 98 ! 99 ENDIF 100 101 iadvtr = iadvtr+1 102 iapptrac = iadvtr 103 104 105 ! Test pour savoir si on advecte a ce pas de temps 106 IF ( iadvtr.EQ.iapp_tracvl ) THEN 107 108 !c .. Modif P.Le Van ( 20/12/97 ) .... 109 !c 110 111 ! traitement des flux de masse avant advection. 112 ! 1. calcul de w 113 ! 2. groupement des mailles pres du pole. 114 115 CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) 116 117 ! ... Flux de masse diaganostiques traceurs 118 flxw = wg / REAL(iapp_tracvl) 119 120 ! test sur l'eventuelle creation de valeurs negatives de la masse 121 DO l=1,llm-1 122 DO ij = iip2+1,ip1jm 123 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & 124 - pbarvg(ij-iip1,l) + pbarvg(ij,l) & 125 + wg(ij,l+1) - wg(ij,l) 126 ENDDO 127 CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 128 DO ij = iip2,ip1jm 129 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 130 ENDDO 131 132 133 CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 134 135 IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN 136 PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, & 137 ' MAX:', zdpmax 138 ENDIF 139 140 ENDDO 141 142 143 !------------------------------------------------------------------- 144 ! Calcul des criteres CFL en X, Y et Z 145 !------------------------------------------------------------------- 146 147 if (countcfl == 0. ) then 148 cflxmax(:)=0. 149 cflymax(:)=0. 150 cflzmax(:)=0. 151 endif 152 153 countcfl=countcfl+iapp_tracvl 154 cflx(:,:)=0. 155 cfly(:,:)=0. 156 cflz(:,:)=0. 157 do l=1,llm 158 do ij=iip2,ip1jm-1 159 if (pbarug(ij,l)>=0.) then 160 cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l) 161 else 162 cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l) 163 endif 164 enddo 165 enddo 166 do l=1,llm 167 do ij=iip2,ip1jm-1,iip1 168 cflx(ij+iip1,l)=cflx(ij,l) 169 enddo 170 enddo 171 172 do l=1,llm 173 do ij=1,ip1jm 174 if (pbarvg(ij,l)>=0.) then 175 cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l) 176 else 177 cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l) 178 endif 179 enddo 180 enddo 181 182 do l=2,llm 183 do ij=1,ip1jm 184 if (wg(ij,l)>=0.) then 185 cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l) 186 else 187 cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1) 188 endif 189 enddo 190 enddo 191 192 do l=1,llm 193 cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l))) 194 cflymax(l)=max(cflymax(l),maxval(cfly(:,l))) 195 cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l))) 196 enddo 203 197 204 198 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 205 ! Par defaut, on sort le diagnostic des CFL tous les jours.206 ! Si on veut le sortir a chaque pas d'advection en cas de plantage207 ! if (countcfl==iapp_tracvl) then199 ! Par defaut, on sort le diagnostic des CFL tous les jours. 200 ! Si on veut le sortir a chaque pas d'advection en cas de plantage 201 ! if (countcfl==iapp_tracvl) then 208 202 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 209 210 211 write(lunout,*) 'L, CFLmax '212 s,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))213 214 215 216 217 c-------------------------------------------------------------------218 cAdvection proprement dite (Modification Le Croller (07/2001)219 c-------------------------------------------------------------------220 221 c----------------------------------------------------222 cCalcul des moyennes basées sur la masse223 c----------------------------------------------------224 225 226 c-----------------------------------------------------------227 cAppel des sous programmes d'advection228 c-----------------------------------------------------------229 230 ccall clock(t_initial)203 if (countcfl==day_step) then 204 do l=1,llm 205 write(lunout,*) 'L, CFLmax ' & 206 ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l)) 207 enddo 208 countcfl=0 209 endif 210 211 !------------------------------------------------------------------- 212 ! Advection proprement dite (Modification Le Croller (07/2001) 213 !------------------------------------------------------------------- 214 215 !---------------------------------------------------- 216 ! Calcul des moyennes basées sur la masse 217 !---------------------------------------------------- 218 call massbar(massem,massebx,masseby) 219 220 !----------------------------------------------------------- 221 ! Appel des sous programmes d'advection 222 !----------------------------------------------------------- 223 do iq=1,nqtot 224 ! call clock(t_initial) 231 225 if(iadv(iq) == 0) cycle 232 c----------------------------------------------------------------233 cSchema de Van Leer I MUSCL234 c----------------------------------------------------------------226 ! ---------------------------------------------------------------- 227 ! Schema de Van Leer I MUSCL 228 ! ---------------------------------------------------------------- 235 229 if(iadv(iq).eq.10) THEN 236 237 238 c----------------------------------------------------------------239 cSchema "pseudo amont" + test sur humidite specifique240 Cpour la vapeur d'eau. F. Codron241 c----------------------------------------------------------------230 call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr) 231 232 ! ---------------------------------------------------------------- 233 ! Schema "pseudo amont" + test sur humidite specifique 234 ! pour la vapeur d'eau. F. Codron 235 ! ---------------------------------------------------------------- 242 236 else if(iadv(iq).eq.14) then 243 c 244 CALL vlspltqs( q(1,1,1), 2., massem, wg , 245 *pbarug,pbarvg,dtvr,p,pk,teta )246 c----------------------------------------------------------------247 cSchema de Frederic Hourdin248 c----------------------------------------------------------------237 ! 238 CALL vlspltqs( q(1,1,1), 2., massem, wg , & 239 pbarug,pbarvg,dtvr,p,pk,teta ) 240 ! ---------------------------------------------------------------- 241 ! Schema de Frederic Hourdin 242 ! ---------------------------------------------------------------- 249 243 else if(iadv(iq).eq.12) then 250 cPas de temps adaptatif244 ! Pas de temps adaptatif 251 245 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 252 246 if (n.GT.1) then 253 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',254 sdtvr,'n=',n247 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 248 dtvr,'n=',n 255 249 endif 256 250 do indice=1,n 257 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)251 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 258 252 end do 259 253 else if(iadv(iq).eq.13) then 260 cPas de temps adaptatif254 ! Pas de temps adaptatif 261 255 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 262 256 if (n.GT.1) then 263 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',264 sdtvr,'n=',n265 endif 266 do indice=1,n267 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)268 end do269 c----------------------------------------------------------------270 cSchema de pente SLOPES271 c----------------------------------------------------------------257 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 258 dtvr,'n=',n 259 endif 260 do indice=1,n 261 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 262 end do 263 ! ---------------------------------------------------------------- 264 ! Schema de pente SLOPES 265 ! ---------------------------------------------------------------- 272 266 else if (iadv(iq).eq.20) then 273 274 275 c----------------------------------------------------------------276 cSchema de Prather277 c----------------------------------------------------------------267 call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 268 269 ! ---------------------------------------------------------------- 270 ! Schema de Prather 271 ! ---------------------------------------------------------------- 278 272 else if (iadv(iq).eq.30) then 279 cPas de temps adaptatif273 ! Pas de temps adaptatif 280 274 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 281 275 if (n.GT.1) then 282 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', 283 s dtvr,'n=',n 284 endif 285 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, 286 s n,dtbon) 287 288 c ---------------------------------------------------------------- 289 c Schemas PPM Lin et Rood 290 c ---------------------------------------------------------------- 291 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. 292 s iadv(iq).LE.18)) then 293 294 c Test sur le flux horizontal 295 c Pas de temps adaptatif 296 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 297 if (n.GT.1) then 298 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', 299 s dtvr,'n=',n 300 endif 301 c Test sur le flux vertical 302 CFLmaxz=0. 303 do l=2,llm 304 do ij=iip2,ip1jm 305 aaa=wg(ij,l)*dtvr/massem(ij,l) 306 CFLmaxz=max(CFLmaxz,aaa) 307 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 308 CFLmaxz=max(CFLmaxz,bbb) 276 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 277 dtvr,'n=',n 278 endif 279 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, & 280 n,dtbon) 281 282 ! ---------------------------------------------------------------- 283 ! Schemas PPM Lin et Rood 284 ! ---------------------------------------------------------------- 285 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. & 286 iadv(iq).LE.18)) then 287 288 ! Test sur le flux horizontal 289 ! Pas de temps adaptatif 290 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 291 if (n.GT.1) then 292 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 293 dtvr,'n=',n 294 endif 295 ! Test sur le flux vertical 296 CFLmaxz=0. 297 do l=2,llm 298 do ij=iip2,ip1jm 299 aaa=wg(ij,l)*dtvr/massem(ij,l) 300 CFLmaxz=max(CFLmaxz,aaa) 301 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 302 CFLmaxz=max(CFLmaxz,bbb) 303 enddo 309 304 enddo 310 enddo 311 if (CFLmaxz.GE.1) then 312 write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 313 endif 314 315 c----------------------------------------------------------- 316 c Ss-prg interface LMDZ.4->PPM3d 317 c----------------------------------------------------------- 318 319 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, 320 s apppm,bpppm,massebx,masseby,pbarug,pbarvg, 321 s unatppm,vnatppm,psppm) 322 323 do indice=1,n 324 c--------------------------------------------------------------------- 325 c VL (version PPM) horiz. et PPM vert. 326 c--------------------------------------------------------------------- 327 if (iadv(iq).eq.11) then 328 c Ss-prg PPM3d de Lin 329 call ppm3d(1,qppm(1,1,iq), 330 s psppm,psppm, 331 s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, 332 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 333 s fill,dum,220.) 334 335 c---------------------------------------------------------------------- 336 c Monotonic PPM 337 c---------------------------------------------------------------------- 338 else if (iadv(iq).eq.16) then 339 c Ss-prg PPM3d de Lin 340 call ppm3d(1,qppm(1,1,iq), 341 s psppm,psppm, 342 s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, 343 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 344 s fill,dum,220.) 345 c--------------------------------------------------------------------- 346 347 c--------------------------------------------------------------------- 348 c Semi Monotonic PPM 349 c--------------------------------------------------------------------- 350 else if (iadv(iq).eq.17) then 351 c Ss-prg PPM3d de Lin 352 call ppm3d(1,qppm(1,1,iq), 353 s psppm,psppm, 354 s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, 355 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 356 s fill,dum,220.) 357 c--------------------------------------------------------------------- 358 359 c--------------------------------------------------------------------- 360 c Positive Definite PPM 361 c--------------------------------------------------------------------- 362 else if (iadv(iq).eq.18) then 363 c Ss-prg PPM3d de Lin 364 call ppm3d(1,qppm(1,1,iq), 365 s psppm,psppm, 366 s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, 367 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 368 s fill,dum,220.) 369 c--------------------------------------------------------------------- 370 endif 371 enddo 372 c----------------------------------------------------------------- 373 c Ss-prg interface PPM3d-LMDZ.4 374 c----------------------------------------------------------------- 375 call interpost(q(1,1,iq),qppm(1,1,iq)) 376 endif 377 c---------------------------------------------------------------------- 378 379 c----------------------------------------------------------------- 380 c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 381 c et Nord j=1 382 c----------------------------------------------------------------- 383 384 c call traceurpole(q(1,1,iq),massem) 385 386 c calcul du temps cpu pour un schema donne 387 388 c call clock(t_final) 389 cym tps_cpu=t_final-t_initial 390 cym cpuadv(iq)=cpuadv(iq)+tps_cpu 391 392 end DO 393 394 395 c------------------------------------------------------------------ 396 c on reinitialise a zero les flux de masse cumules 397 c--------------------------------------------------- 398 iadvtr=0 399 400 ENDIF ! if iadvtr.EQ.iapp_tracvl 401 402 RETURN 403 END 404 305 if (CFLmaxz.GE.1) then 306 write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 307 endif 308 309 !----------------------------------------------------------- 310 ! Ss-prg interface LMDZ.4->PPM3d 311 !----------------------------------------------------------- 312 313 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & 314 apppm,bpppm,massebx,masseby,pbarug,pbarvg, & 315 unatppm,vnatppm,psppm) 316 317 do indice=1,n 318 !---------------------------------------------------------------- 319 ! VL (version PPM) horiz. et PPM vert. 320 !---------------------------------------------------------------- 321 if (iadv(iq).eq.11) then 322 ! Ss-prg PPM3d de Lin 323 call ppm3d(1,qppm(1,1,iq), & 324 psppm,psppm, & 325 unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, & 326 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 327 fill,dum,220.) 328 329 !------------------------------------------------------------- 330 ! Monotonic PPM 331 !------------------------------------------------------------- 332 else if (iadv(iq).eq.16) then 333 ! Ss-prg PPM3d de Lin 334 call ppm3d(1,qppm(1,1,iq), & 335 psppm,psppm, & 336 unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, & 337 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 338 fill,dum,220.) 339 !------------------------------------------------------------- 340 341 !------------------------------------------------------------- 342 ! Semi Monotonic PPM 343 !------------------------------------------------------------- 344 else if (iadv(iq).eq.17) then 345 ! Ss-prg PPM3d de Lin 346 call ppm3d(1,qppm(1,1,iq), & 347 psppm,psppm, & 348 unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, & 349 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 350 fill,dum,220.) 351 !------------------------------------------------------------- 352 353 !------------------------------------------------------------- 354 ! Positive Definite PPM 355 !------------------------------------------------------------- 356 else if (iadv(iq).eq.18) then 357 ! Ss-prg PPM3d de Lin 358 call ppm3d(1,qppm(1,1,iq), & 359 psppm,psppm, & 360 unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, & 361 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 362 fill,dum,220.) 363 !------------------------------------------------------------- 364 endif 365 enddo 366 !----------------------------------------------------------------- 367 ! Ss-prg interface PPM3d-LMDZ.4 368 !----------------------------------------------------------------- 369 call interpost(q(1,1,iq),qppm(1,1,iq)) 370 endif 371 !---------------------------------------------------------------------- 372 373 !----------------------------------------------------------------- 374 ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 375 ! et Nord j=1 376 !----------------------------------------------------------------- 377 378 ! call traceurpole(q(1,1,iq),massem) 379 380 ! calcul du temps cpu pour un schema donne 381 382 ! call clock(t_final) 383 !ym tps_cpu=t_final-t_initial 384 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu 385 386 end DO 387 388 389 !------------------------------------------------------------------ 390 ! on reinitialise a zero les flux de masse cumules 391 !--------------------------------------------------- 392 iadvtr=0 393 394 ENDIF ! if iadvtr.EQ.iapp_tracvl 395 396 END SUBROUTINE advtrac -
trunk/LMDZ.COMMON/libf/dyn3d/ce0l.F90
r101 r270 1 1 ! 2 ! $Id: ce0l.F90 1 425 2010-09-02 13:45:23Z lguez$2 ! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $ 3 3 ! 4 4 !------------------------------------------------------------------------------- … … 27 27 #endif 28 28 IMPLICIT NONE 29 #include "iniprint.h"30 29 #ifndef CPP_EARTH 31 30 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' … … 37 36 #include "paramet.h" 38 37 #include "indicesol.h" 38 #include "iniprint.h" 39 39 #include "temps.h" 40 40 #include "logic.h" 41 INTEGER, PARAMETER :: longcles=20 42 REAL, DIMENSION(longcles) :: clesphy0 41 43 REAL, DIMENSION(iip1,jjp1) :: masque 42 44 CHARACTER(LEN=15) :: calnd 45 REAL, DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol 43 46 !------------------------------------------------------------------------------- 44 CALL conf_gcm( 99, .TRUE. )47 CALL conf_gcm( 99, .TRUE. , clesphy0 ) 45 48 46 49 CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) … … 78 81 WRITE(lunout,'(//)') 79 82 WRITE(lunout,*) ' interbar = ',interbar 80 CALL etat0_netcdf(interbar,masque, ok_etat0)83 CALL etat0_netcdf(interbar,masque,phis,ok_etat0) 81 84 82 85 IF(ok_limit) THEN … … 89 92 END IF 90 93 94 IF (grilles_gcm_netcdf) THEN 95 WRITE(lunout,'(//)') 96 WRITE(lunout,*) ' *************************** ' 97 WRITE(lunout,*) ' *** grilles_gcm_netcdf *** ' 98 WRITE(lunout,*) ' *************************** ' 99 WRITE(lunout,'(//)') 100 CALL grilles_gcm_netcdf_sub(masque,phis) 101 END IF 91 102 #endif 92 103 ! of #ifndef CPP_EARTH #else -
trunk/LMDZ.COMMON/libf/dyn3d/comdissipn.h
r1 r270 2 2 ! $Header$ 3 3 ! 4 c----------------------------------------------------------------------- 5 c INCLUDE comdissipn.h 4 ! Attention : ce fichier include est compatible format fixe/format libre 5 ! veillez à n'utiliser que des ! pour les commentaires 6 ! et à bien positionner les & des lignes de continuation 7 ! (les placer en colonne 6 et en colonne 73) 8 !----------------------------------------------------------------------- 9 ! INCLUDE comdissipn.h 6 10 7 11 REAL tetaudiv, tetaurot, tetah, cdivu, crot, cdivh 8 c 9 COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm) , 10 1cdivu, crot, cdivh12 ! 13 COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm) , & 14 & cdivu, crot, cdivh 11 15 12 c 13 cLes parametres de ce common proviennent des calculs effectues dans14 cInidissip .15 c 16 c-----------------------------------------------------------------------16 ! 17 ! Les parametres de ce common proviennent des calculs effectues dans 18 ! Inidissip . 19 ! 20 !----------------------------------------------------------------------- -
trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F
r119 r270 230 230 CALL getin('output_grads_dyn',output_grads_dyn) 231 231 232 !Config Key = idissip232 !Config Key = dissip_period 233 233 !Config Desc = periode de la dissipation 234 !Config Def = 10234 !Config Def = 0 235 235 !Config Help = periode de la dissipation 236 !Config (en pas) ... a completer ! 237 idissip = 10 238 CALL getin('idissip',idissip) 236 !Config dissip_period=0 => la valeur sera calcule dans inidissip 237 !Config dissip_period>0 => on prend cette valeur 238 dissip_period = 0 239 CALL getin('dissip_period',dissip_period) 239 240 240 241 ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... … … 652 653 write(lunout,*)' periodav = ', periodav 653 654 write(lunout,*)' output_grads_dyn = ', output_grads_dyn 654 write(lunout,*)' idissip = ', idissip655 write(lunout,*)' dissip_period = ', dissip_period 655 656 write(lunout,*)' lstardis = ', lstardis 656 657 write(lunout,*)' nitergdiv = ', nitergdiv … … 876 877 CALL getin('ok_etat0',ok_etat0) 877 878 879 !Config Key = grilles_gcm_netcdf 880 !Config Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit 881 !Config Def = n 882 grilles_gcm_netcdf = .FALSE. 883 CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf) 884 878 885 write(lunout,*)' #########################################' 879 886 write(lunout,*)' Configuration des parametres de cel0' … … 890 897 write(lunout,*)' periodav = ', periodav 891 898 write(lunout,*)' output_grads_dyn = ', output_grads_dyn 892 write(lunout,*)' idissip = ', idissip899 write(lunout,*)' dissip_period = ', dissip_period 893 900 write(lunout,*)' lstardis = ', lstardis 894 901 write(lunout,*)' nitergdiv = ', nitergdiv … … 922 929 write(lunout,*)' ok_limit = ', ok_limit 923 930 write(lunout,*)' ok_etat0 = ', ok_etat0 931 write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf 924 932 c 925 933 RETURN -
trunk/LMDZ.COMMON/libf/dyn3d/control_mod.F90
r97 r270 12 12 REAL :: periodav 13 13 INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys 14 INTEGER :: iconser,iecri, idissip,iphysiq,iecrimoy14 INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy 15 15 INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn 16 16 LOGICAL :: offline -
trunk/LMDZ.COMMON/libf/dyn3d/defrun.F
r6 r270 139 139 140 140 READ (tapedef,9001) ch1,ch4 141 READ (tapedef,*) idissip142 WRITE(tapeout,9001) ch1,' idissip'143 WRITE(tapeout,*) idissip141 READ (tapedef,*) dissip_period 142 WRITE(tapeout,9001) ch1,'dissip_period' 143 WRITE(tapeout,*) dissip_period 144 144 145 145 ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... -
trunk/LMDZ.COMMON/libf/dyn3d/gcm.F
r130 r270 427 427 c ------------------------------- 428 428 429 IF (call_iniphys.and.(iflag_phys .eq.1)) THEN429 IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN 430 430 latfi(1)=rlatu(1) 431 431 lonfi(1)=0. … … 486 486 #endif 487 487 488 #ifdef CPP_EARTH 489 ! Create start file (startphy.nc) and boundary conditions (limit.nc) 490 ! for the Earth verstion 491 if (iflag_phys>=100) then 492 call iniaqua(ngridmx,latfi,lonfi,iflag_phys) 493 endif 494 #endif 495 488 496 if (planet_type.eq."mars") then 489 497 ! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars -
trunk/LMDZ.COMMON/libf/dyn3d/ini_paramLMDZ_dyn.h
r1 r270 78 78 . "once", dt_cum,dt_cum) 79 79 c 80 CALL histdef(nid_ctesGCM, " idissip",80 CALL histdef(nid_ctesGCM, "dissip_period", 81 81 . "periode de la dissipation (en pas) ... a completer", 82 82 . "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, -
trunk/LMDZ.COMMON/libf/dyn3d/iniacademic.F90
r127 r270 1 1 ! 2 ! $Id: iniacademic.F90 152 0 2011-05-23 11:37:09Z emillour$2 ! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $ 3 3 ! 4 4 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) … … 115 115 endif 116 116 117 academic_case: if (iflag_phys == 2) then117 academic_case: if (iflag_phys >= 2) then 118 118 ! initializations 119 119 … … 208 208 IF (.NOT. read_start) THEN 209 209 ! surface pressure 210 ps(:)=preff 210 if (iflag_phys>2) then 211 ps(:)=108080. ! Earth aqua/terra planets 212 else 213 ps(:)=preff 214 endif 211 215 ! ground geopotential 212 216 phis(:)=0. -
trunk/LMDZ.COMMON/libf/dyn3d/iniconst.F
r127 r270 47 47 c----------------------------------------------------------------------- 48 48 49 dtdiss = idissip * dtvr50 49 dtphys = iphysiq * dtvr 51 50 unsim = 1./iim -
trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90
r269 r270 1 1 ! 2 ! $Id: inidissip.F 1403 2010-07-01 09:02:53Z fairhead$2 ! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $ 3 3 ! 4 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , 5 * tetagdiv,tetagrot,tetatemp ) 6 c======================================================================= 7 c initialisation de la dissipation horizontale 8 c======================================================================= 9 c----------------------------------------------------------------------- 10 c declarations: 11 c ------------- 12 13 USE control_mod 14 15 IMPLICIT NONE 16 #include "dimensions.h" 17 #include "paramet.h" 18 #include "comdissipn.h" 19 #include "comconst.h" 20 #include "comvert.h" 21 #include "logic.h" 22 23 LOGICAL lstardis 24 INTEGER nitergdiv,nitergrot,niterh 25 REAL tetagdiv,tetagrot,tetatemp 26 REAL fact,zvert(llm),zz 27 REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm) 28 REAL ullm,vllm,umin,vmin,zhmin,zhmax 29 REAL zllm,z1llm 30 31 INTEGER l,ij,idum,ii 32 REAL tetamin 33 REAL Pup 34 35 REAL ran1 36 37 38 c----------------------------------------------------------------------- 39 c 40 c calcul des valeurs propres des operateurs par methode iterrative: 41 c ----------------------------------------------------------------- 42 43 crot = -1. 44 cdivu = -1. 45 cdivh = -1. 46 47 c calcul de la valeur propre de divgrad: 48 c -------------------------------------- 49 idum = 0 50 DO l = 1, llm 51 DO ij = 1, ip1jmp1 4 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , & 5 tetagdiv,tetagrot,tetatemp ) 6 !======================================================================= 7 ! initialisation de la dissipation horizontale 8 !======================================================================= 9 !----------------------------------------------------------------------- 10 ! declarations: 11 ! ------------- 12 13 USE control_mod, only : dissip_period,iperiod 14 15 IMPLICIT NONE 16 include "dimensions.h" 17 include "paramet.h" 18 include "comdissipn.h" 19 include "comconst.h" 20 include "comvert.h" 21 include "logic.h" 22 include "iniprint.h" 23 24 LOGICAL,INTENT(in) :: lstardis 25 INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh 26 REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp 27 28 ! Local variables: 29 REAL fact,zvert(llm),zz 30 REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm) 31 REAL ullm,vllm,umin,vmin,zhmin,zhmax 32 REAL zllm,z1llm 33 34 INTEGER l,ij,idum,ii 35 REAL tetamin 36 REAL Pup 37 character (len=80) :: abort_message 38 39 REAL ran1 40 41 42 !----------------------------------------------------------------------- 43 ! 44 ! calcul des valeurs propres des operateurs par methode iterrative: 45 ! ----------------------------------------------------------------- 46 47 crot = -1. 48 cdivu = -1. 49 cdivh = -1. 50 51 ! calcul de la valeur propre de divgrad: 52 ! -------------------------------------- 53 idum = 0 54 DO l = 1, llm 55 DO ij = 1, ip1jmp1 52 56 deltap(ij,l) = 1. 53 ENDDO 54 ENDDO 55 56 idum = -1 57 zh(1) = RAN1(idum)-.5 58 idum = 0 59 DO ij = 2, ip1jmp1 60 zh(ij) = RAN1(idum) -.5 61 ENDDO 62 63 CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) 64 65 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 66 67 IF ( zhmin .GE. zhmax ) THEN 68 PRINT*,' Inidissip zh min max ',zhmin,zhmax 69 STOP'probleme generateur alleatoire dans inidissip' 70 ENDIF 71 72 zllm = ABS( zhmax ) 73 DO l = 1,50 74 IF(lstardis) THEN 75 CALL divgrad2(1,zh,deltap,niterh,zh) 76 ELSE 77 CALL divgrad (1,zh,niterh,zh) 78 ENDIF 79 80 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 81 82 zllm = ABS( zhmax ) 83 z1llm = 1./zllm 84 DO ij = 1,ip1jmp1 85 zh(ij) = zh(ij)* z1llm 86 ENDDO 87 ENDDO 88 89 IF(lstardis) THEN 90 cdivh = 1./ zllm 91 ELSE 92 cdivh = zllm ** ( -1./niterh ) 93 ENDIF 94 95 c calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 96 c ----------------------------------------------------------------- 97 print*,'calcul des valeurs propres' 98 99 DO 20 ii = 1, 2 100 c 101 DO ij = 1, ip1jmp1 102 zu(ij) = RAN1(idum) -.5 103 ENDDO 104 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 105 DO ij = 1, ip1jm 106 zv(ij) = RAN1(idum) -.5 107 ENDDO 108 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 109 110 CALL minmax(iip1*jjp1,zu,umin,ullm ) 111 CALL minmax(iip1*jjm, zv,vmin,vllm ) 112 113 ullm = ABS ( ullm ) 114 vllm = ABS ( vllm ) 115 116 DO 5 l = 1, 50 117 IF(ii.EQ.1) THEN 118 ccccc CALL covcont( 1,zu,zv,zu,zv ) 119 IF(lstardis) THEN 120 CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv ) 121 ELSE 122 CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv ) 123 ENDIF 124 ELSE 125 IF(lstardis) THEN 126 CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv ) 127 ELSE 128 CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv ) 129 ENDIF 130 ENDIF 131 132 CALL minmax(iip1*jjp1,zu,umin,ullm ) 133 CALL minmax(iip1*jjm, zv,vmin,vllm ) 134 135 ullm = ABS ( ullm ) 136 vllm = ABS ( vllm ) 137 138 zllm = MAX( ullm,vllm ) 139 z1llm = 1./ zllm 140 DO ij = 1, ip1jmp1 141 zu(ij) = zu(ij)* z1llm 142 ENDDO 143 DO ij = 1, ip1jm 144 zv(ij) = zv(ij)* z1llm 145 ENDDO 146 5 CONTINUE 147 148 IF ( ii.EQ.1 ) THEN 149 IF(lstardis) THEN 150 cdivu = 1./zllm 151 ELSE 152 cdivu = zllm **( -1./nitergdiv ) 153 ENDIF 154 ELSE 155 IF(lstardis) THEN 156 crot = 1./ zllm 157 ELSE 158 crot = zllm **( -1./nitergrot ) 159 ENDIF 160 ENDIF 161 162 20 CONTINUE 163 164 c petit test pour les operateurs non star: 165 c ---------------------------------------- 166 167 c IF(.NOT.lstardis) THEN 168 fact = rad*24./REAL(jjm) 169 fact = fact*fact 170 PRINT*,'coef u ', fact/cdivu, 1./cdivu 171 PRINT*,'coef r ', fact/crot , 1./crot 172 PRINT*,'coef h ', fact/cdivh, 1./cdivh 173 c ENDIF 174 175 c----------------------------------------------------------------------- 176 c variation verticale du coefficient de dissipation: 177 c -------------------------------------------------- 178 179 c First step: going from 1 to dissip_fac_mid (in gcm.def) 180 c============ 181 DO l=1,llm 182 zz = 1. - preff/presnivs(l) 183 zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) 184 ENDDO 185 186 write(*,*) 'Dissipation : ' 187 write(*,*) 'Multiplication de la dissipation en altitude :' 188 write(*,*) ' dissip_fac_mid =', dissip_fac_mid 189 190 c Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) 191 c========================== 192 c Utilisation de la fonction tangente hyperbolique pour augmenter 193 c arbitrairement la dissipation et donc la stabilite du modele a 194 c partir d'une certaine altitude. 195 196 c Le facteur multiplicatif de basse atmosphere etant deja pris 197 c en compte, il faut diviser le facteur multiplicatif de haute 198 c atmosphere par celui-ci. 199 200 if (ok_strato) then 201 202 Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) 203 do l=1,llm 204 zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) 205 & *(1-(0.5*(1+tanh(-6./dissip_deltaz* 206 & (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) 207 enddo 208 209 write(*,*) ' dissip_fac_up =', dissip_fac_up 210 write(*,*) 'Transition mid /up: Pupstart,delta =', 211 & dissip_pupstart,'Pa', dissip_deltaz , '(km)' 212 213 endif 214 215 216 PRINT*,'Constantes de temps de la diffusion horizontale' 217 218 tetamin = 1.e+6 219 220 DO l=1,llm 221 tetaudiv(l) = zvert(l)/tetagdiv 222 tetaurot(l) = zvert(l)/tetagrot 223 tetah(l) = zvert(l)/tetatemp 224 225 IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 226 IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 227 IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) 228 ENDDO 229 230 PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod 231 idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod 232 PRINT *,' INIDI tetamin idissip ',tetamin,idissip 233 idissip = MAX(iperiod,idissip) 234 dtdiss = idissip * dtvr 235 PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss 236 237 DO l = 1,llm 238 PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), 239 * dtdiss*tetah(l) 240 ENDDO 241 242 c 243 RETURN 244 END 57 ENDDO 58 ENDDO 59 60 idum = -1 61 zh(1) = RAN1(idum)-.5 62 idum = 0 63 DO ij = 2, ip1jmp1 64 zh(ij) = RAN1(idum) -.5 65 ENDDO 66 67 CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) 68 69 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 70 71 IF ( zhmin .GE. zhmax ) THEN 72 write(lunout,*)' Inidissip zh min max ',zhmin,zhmax 73 abort_message='probleme generateur alleatoire dans inidissip' 74 call abort_gcm('inidissip',abort_message,1) 75 ENDIF 76 77 zllm = ABS( zhmax ) 78 DO l = 1,50 79 IF(lstardis) THEN 80 CALL divgrad2(1,zh,deltap,niterh,zh) 81 ELSE 82 CALL divgrad (1,zh,niterh,zh) 83 ENDIF 84 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 92 ENDDO 93 94 IF(lstardis) THEN 95 cdivh = 1./ zllm 96 ELSE 97 cdivh = zllm ** ( -1./niterh ) 98 ENDIF 99 100 ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 101 ! ----------------------------------------------------------------- 102 write(lunout,*)'inidissip: calcul des valeurs propres' 103 104 DO ii = 1, 2 105 ! 106 DO ij = 1, ip1jmp1 107 zu(ij) = RAN1(idum) -.5 108 ENDDO 109 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 110 DO ij = 1, ip1jm 111 zv(ij) = RAN1(idum) -.5 112 ENDDO 113 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 114 115 CALL minmax(iip1*jjp1,zu,umin,ullm ) 116 CALL minmax(iip1*jjm, zv,vmin,vllm ) 117 118 ullm = ABS ( ullm ) 119 vllm = ABS ( vllm ) 120 121 DO l = 1, 50 122 IF(ii.EQ.1) THEN 123 !cccc CALL covcont( 1,zu,zv,zu,zv ) 124 IF(lstardis) THEN 125 CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv ) 126 ELSE 127 CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv ) 128 ENDIF 129 ELSE 130 IF(lstardis) THEN 131 CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv ) 132 ELSE 133 CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv ) 134 ENDIF 135 ENDIF 136 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 151 end DO 152 153 IF ( ii.EQ.1 ) THEN 154 IF(lstardis) THEN 155 cdivu = 1./zllm 156 ELSE 157 cdivu = zllm **( -1./nitergdiv ) 158 ENDIF 159 ELSE 160 IF(lstardis) THEN 161 crot = 1./ zllm 162 ELSE 163 crot = zllm **( -1./nitergrot ) 164 ENDIF 165 ENDIF 166 167 end DO 168 169 ! petit test pour les operateurs non star: 170 ! ---------------------------------------- 171 172 ! IF(.NOT.lstardis) THEN 173 fact = rad*24./REAL(jjm) 174 fact = fact*fact 175 write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu 176 write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot 177 write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh 178 ! ENDIF 179 180 !----------------------------------------------------------------------- 181 ! variation verticale du coefficient de dissipation: 182 ! -------------------------------------------------- 183 184 ! First step: going from 1 to dissip_fac_mid (in gcm.def) 185 !============ 186 DO l=1,llm 187 zz = 1. - preff/presnivs(l) 188 zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) 189 ENDDO 190 191 write(lunout,*) 'Dissipation : ' 192 write(lunout,*) 'Multiplication de la dissipation en altitude :' 193 write(lunout,*) ' dissip_fac_mid =', dissip_fac_mid 194 195 ! Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) 196 !========================== 197 ! Utilisation de la fonction tangente hyperbolique pour augmenter 198 ! arbitrairement la dissipation et donc la stabilite du modele a 199 ! partir d'une certaine altitude. 200 201 ! Le facteur multiplicatif de basse atmosphere etant deja pris 202 ! en compte, il faut diviser le facteur multiplicatif de haute 203 ! atmosphere par celui-ci. 204 205 if (ok_strato) then 206 207 Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) 208 do l=1,llm 209 zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) & 210 *(1-(0.5*(1+tanh(-6./dissip_deltaz* & 211 (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) 212 enddo 213 214 write(*,*) ' dissip_fac_up =', dissip_fac_up 215 write(*,*) 'Transition mid /up: Pupstart,delta =', & 216 dissip_pupstart,'Pa', dissip_deltaz , '(km)' 217 218 endif 219 220 221 write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' 222 223 tetamin = 1.e+6 224 225 DO l=1,llm 226 tetaudiv(l) = zvert(l)/tetagdiv 227 tetaurot(l) = zvert(l)/tetagrot 228 tetah(l) = zvert(l)/tetatemp 229 230 IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 231 IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 232 IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) 233 ENDDO 234 235 ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def 236 IF (dissip_period == 0) THEN 237 dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod 238 write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period 239 dissip_period = MAX(iperiod,dissip_period) 240 END IF 241 242 dtdiss = dissip_period * dtvr 243 write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr 244 245 DO l = 1,llm 246 write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & 247 dtdiss*tetah(l) 248 ENDDO 249 250 END SUBROUTINE inidissip -
trunk/LMDZ.COMMON/libf/dyn3d/integrd.F
r7 r270 1 1 ! 2 ! $Id: integrd.F 1 446 2010-10-22 09:27:25Z emillour$2 ! $Id: integrd.F 1550 2011-07-05 09:44:55Z lguez $ 3 3 ! 4 4 SUBROUTINE integrd … … 89 89 IF( ps(ij).LT.0. ) THEN 90 90 PRINT*,' Au point ij = ',ij, ' , pression sol neg. ', ps(ij) 91 STOP' dans integrd' 91 print *, ' dans integrd' 92 stop 1 92 93 ENDIF 93 94 ENDDO -
trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F
r130 r270 152 152 logical ok_sync 153 153 parameter (ok_sync = .true.) 154 logical physic 154 155 155 156 data callinigrads/.true./ … … 222 223 223 224 itau = 0 225 physic=.true. 226 if (iflag_phys==0.or.iflag_phys==2) physic=.false. 227 224 228 c iday = day_ini+itau/day_step 225 229 c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0 … … 322 326 ! Purely Matsuno time stepping 323 327 IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE. 324 IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE. 328 IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward ) 329 s apdiss = .TRUE. 325 330 IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 326 s .and. iflag_phys.EQ.1) apphys = .TRUE.331 s .and. physic ) apphys = .TRUE. 327 332 ELSE 328 333 ! Leapfrog/Matsuno time stepping 329 334 IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE. 330 IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE. 331 IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE. 335 IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward ) 336 s apdiss = .TRUE. 337 IF( MOD(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.TRUE. 332 338 END IF 333 339 … … 357 363 c ------------------------------------------------------------- 358 364 359 IF( forward. OR . leapf ) THEN 360 365 ! IF( forward. OR . leapf ) THEN 366 IF((.not.forward).OR. leapf ) THEN 367 ! Ehouarn: gather mass fluxes during backward Matsuno or LF step 361 368 CALL caladvtrac(q,pbaru,pbarv, 362 369 * p, masse, dq, teta, … … 497 504 ENDDO 498 505 ENDDO ! of DO l=1,llm 499 call friction(ucov,vcov,dtvr)500 506 501 if (planet_type.eq."giant") then 502 ! Intrinsic heat flux 503 ! Aymeric -- for giant planets 504 if (ihf .gt. 1.e-6) then 505 !print *, '**** INTRINSIC HEAT FLUX ****', ihf 506 DO ij=1,ip1jmp1 507 teta(ij,1) = teta(ij,1) 508 & + dtvr * aire(ij) * ihf / cpp / masse(ij,1) 509 ENDDO 510 !print *, '**** d teta ' 511 !print *, dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1) 512 endif 513 endif 507 if (planet_type.eq."giant") then 508 ! add an intrinsic heat flux at the base of the atmosphere 509 teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1) 510 endif 511 512 call friction(ucov,vcov,dtvr) 514 513 515 514 -
trunk/LMDZ.COMMON/libf/dyn3d/limit_netcdf.F90
r7 r270 1 1 ! 2 ! $Id: limit_netcdf.F90 1 441 2010-10-13 13:06:56Z emillour$2 ! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $ 3 3 !------------------------------------------------------------------------------- 4 4 ! … … 42 42 REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask 43 43 #ifndef CPP_EARTH 44 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'44 CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1) 45 45 #else 46 46 !------------------------------------------------------------------------------- … … 52 52 #include "indicesol.h" 53 53 54 !--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------55 LOGICAL, PARAMETER :: fracterre=.TRUE.56 57 54 !--- INPUT NETCDF FILES NAMES -------------------------------------------------- 58 55 CHARACTER(LEN=25) :: icefile, sstfile, dumstr 59 56 CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc ', & 60 57 famipsic='amipbc_sic_1x1.nc ', & 61 fclimsst='amipbc_sst_1x1_clim.nc ', &62 fclimsic='amipbc_sic_1x1_clim.nc ', &63 58 fcpldsst='cpl_atm_sst.nc ', & 64 59 fcpldsic='cpl_atm_sic.nc ', & 60 fhistsst='histmth_sst.nc ', & 61 fhistsic='histmth_sic.nc ', & 65 62 frugo ='Rugos.nc ', & 66 63 falbe ='Albedo.nc ' 67 64 CHARACTER(LEN=10) :: varname 68 65 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------ 69 66 REAL, DIMENSION(klon) :: fi_ice, verif … … 80 77 INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC 81 78 INTEGER :: NF90_FORMAT 82 LOGICAL :: lCPL !--- T: IPCC-IPSL cpl model output files83 79 INTEGER :: ndays !--- Depending on the output calendar 84 80 … … 104 100 105 101 !--- RUGOSITY TREATMENT -------------------------------------------------------- 106 WRITE(lunout,*) 'Traitement de la rugosite' 107 CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 102 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite' 103 varname='RUGOS' 104 CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 108 105 109 106 !--- OCEAN TREATMENT ----------------------------------------------------------- 110 PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.107 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique' 111 108 112 109 ! Input SIC file selection 113 icefile='fake' 114 IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic) 115 IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic) 116 IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic) 117 SELECT CASE(icefile) 118 CASE(famipsic); dumstr='Amip.' 119 CASE(fclimsic); dumstr='Amip climatologique.' 120 CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. 121 CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1) 122 END SELECT 110 ! Open file only to test if available 111 IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 112 icefile=TRIM(famipsic) 113 varname='sicbcs' 114 ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 115 icefile=TRIM(fcpldsic) 116 varname='SIICECOV' 117 ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 118 icefile=TRIM(fhistsic) 119 varname='pourc_sic' 120 ELSE 121 WRITE(lunout,*) 'ERROR! No sea-ice input file was found.' 122 WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic) 123 CALL abort_gcm('limit_netcdf','No sea-ice file was found',1) 124 END IF 123 125 ierr=NF90_CLOSE(nid) 124 WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)125 126 CALL get_2Dfield(icefile, 'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL)126 IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile) 127 128 CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice) 127 129 128 130 ALLOCATE(pctsrf_t(klon,nbsrf,ndays)) 129 131 DO k=1,ndays 130 fi_ice=phy_ice(:,k) 131 WHERE(fi_ice>=1.0 ) fi_ice=1.0 132 WHERE(fi_ice<EPSFRA) fi_ice=0.0 133 IF(fracterre) THEN 134 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil 135 pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice 136 IF(lCPL) THEN ! SIC=pICE*(1-LIC-TER) 137 pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 138 ELSE ! SIC=pICE-LIC 132 fi_ice=phy_ice(:,k) 133 WHERE(fi_ice>=1.0 ) fi_ice=1.0 134 WHERE(fi_ice<EPSFRA) fi_ice=0.0 135 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil 136 pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice 137 IF (icefile==trim(fcpldsic)) THEN ! SIC=pICE*(1-LIC-TER) 138 pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 139 ELSE IF (icefile==trim(fhistsic)) THEN ! SIC=pICE 140 pctsrf_t(:,is_sic,k)=fi_ice(:) 141 ELSE ! icefile==famipsic ! SIC=pICE-LIC 139 142 pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k) 140 141 142 143 END IF 144 WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0. 145 WHERE(1.0-zmasq<EPSFRA) 143 146 pctsrf_t(:,is_sic,k)=0.0 144 147 pctsrf_t(:,is_oce,k)=0.0 145 148 ELSEWHERE 146 149 WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq) 147 pctsrf_t(:,is_sic,k)=1.0-zmasq148 pctsrf_t(:,is_oce,k)=0.0150 pctsrf_t(:,is_sic,k)=1.0-zmasq 151 pctsrf_t(:,is_oce,k)=0.0 149 152 ELSEWHERE 150 pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)151 WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)152 pctsrf_t(:,is_oce,k)=0.0153 pctsrf_t(:,is_sic,k)=1.0-zmasq154 END WHERE153 pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k) 154 WHERE(pctsrf_t(:,is_oce,k)<EPSFRA) 155 pctsrf_t(:,is_oce,k)=0.0 156 pctsrf_t(:,is_sic,k)=1.0-zmasq 157 END WHERE 155 158 END WHERE 156 END WHERE 157 nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) 158 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 159 nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) 160 IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad 161 ELSE 162 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) 163 WHERE(NINT(pctsrf(:,is_ter))==1) 164 pctsrf_t(:,is_sic,k)=0. 165 pctsrf_t(:,is_oce,k)=0. 166 WHERE(fi_ice>=1.e-5) 167 pctsrf_t(:,is_lic,k)=fi_ice 168 pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k) 169 ELSEWHERE 170 pctsrf_t(:,is_lic,k)=0.0 171 END WHERE 172 ELSEWHERE 173 pctsrf_t(:,is_lic,k) = 0.0 174 WHERE(fi_ice>=1.e-5) 175 pctsrf_t(:,is_ter,k)=0.0 176 pctsrf_t(:,is_sic,k)=fi_ice 177 pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k) 178 ELSEWHERE 179 pctsrf_t(:,is_sic,k)=0.0 180 pctsrf_t(:,is_oce,k)=1.0 181 END WHERE 182 END WHERE 183 verif=sum(pctsrf_t(:,:,k),dim=2) 184 nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5) 185 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 186 END IF 159 END WHERE 160 nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) 161 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 162 nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) 163 IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad 187 164 END DO 188 165 DEALLOCATE(phy_ice) 189 166 190 167 !--- SST TREATMENT ------------------------------------------------------------- 191 WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.168 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst' 192 169 193 170 ! Input SST file selection 194 sstfile='fake' 195 IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst) 196 IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst) 197 IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst) 198 SELECT CASE(icefile) 199 CASE(famipsic); dumstr='Amip.' 200 CASE(fclimsic); dumstr='Amip climatologique.' 201 CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. 202 CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1) 203 END SELECT 171 ! Open file only to test if available 172 IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 173 sstfile=TRIM(famipsst) 174 varname='tosbcs' 175 ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 176 sstfile=TRIM(fcpldsst) 177 varname='SISUTESW' 178 ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 179 sstfile=TRIM(fhistsst) 180 varname='tsol_oce' 181 ELSE 182 WRITE(lunout,*) 'ERROR! No sst input file was found.' 183 WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst) 184 CALL abort_gcm('limit_netcdf','No sst file was found',1) 185 END IF 204 186 ierr=NF90_CLOSE(nid) 205 WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)206 207 CALL get_2Dfield( trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)187 IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile) 188 189 CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap) 208 190 209 191 !--- ALBEDO TREATMENT ---------------------------------------------------------- 210 WRITE(lunout,*) 'Traitement de l albedo' 211 CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb) 192 IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo' 193 varname='ALBEDO' 194 CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb) 212 195 213 196 !--- REFERENCE GROUND HEAT FLUX TREATMENT -------------------------------------- … … 215 198 216 199 !--- OUTPUT FILE WRITING ------------------------------------------------------- 217 WRITE(lunout,*) 'Ecriture du fichier limit'200 IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut' 218 201 219 202 !--- File creation … … 264 247 ierr=NF90_CLOSE(nid) 265 248 249 IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin' 250 266 251 DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug) 267 252 … … 276 261 !------------------------------------------------------------------------------- 277 262 ! 278 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)263 SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask) 279 264 ! 280 265 !----------------------------------------------------------------------------- … … 304 289 ! Arguments: 305 290 CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name 291 CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name 306 292 CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB 307 293 LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels 308 294 INTEGER, INTENT(IN) :: ndays ! current year number of days 309 REAL, POINTER, DIMENSION(:, :) :: champo 295 REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) 310 296 LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC) 311 297 REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask 312 LOGICAL, OPTIONAL, INTENT(IN) :: lCPL ! Coupled model flag (for ICE)313 298 !------------------------------------------------------------------------------ 314 299 ! Local variables: … … 316 301 INTEGER :: ncid, varid ! NetCDF identifiers 317 302 CHARACTER(LEN=30) :: dnam ! dimension name 318 CHARACTER(LEN=80) :: varname ! NetCDF variable name319 303 !--- dimensions 320 304 INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers … … 331 315 !--- input files 332 316 CHARACTER(LEN=20) :: cal_in ! calendar 317 CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file 333 318 INTEGER :: ndays_in ! number of days 334 319 !--- misc … … 337 322 CHARACTER(LEN=25) :: title ! for messages 338 323 LOGICAL :: extrp ! flag for extrapolation 324 LOGICAL :: oldice ! flag for old way ice computation 339 325 REAL :: chmin, chmax 340 326 INTEGER ierr 341 327 integer n_extrap ! number of extrapolated points 342 328 logical skip 329 343 330 !------------------------------------------------------------------------------ 344 331 !---Variables depending on keyword 'mode' ------------------------------------- 345 332 NULLIFY(champo) 333 346 334 SELECT CASE(mode) 347 CASE('RUG'); varname='RUGOS';title='Rugosite'348 CASE('SIC'); varname='sicbcs'; title='Sea-ice'349 CASE('SST'); varname='tosbcs'; title='SST'350 CASE('ALB'); varname='ALBEDO'; title='Albedo'335 CASE('RUG'); title='Rugosite' 336 CASE('SIC'); title='Sea-ice' 337 CASE('SST'); title='SST' 338 CASE('ALB'); title='Albedo' 351 339 END SELECT 340 341 352 342 extrp=.FALSE. 343 oldice=.FALSE. 353 344 IF ( PRESENT(flag) ) THEN 354 345 IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 346 IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 355 347 END IF 356 348 357 349 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- 350 IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam 358 351 ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid); CALL ncerr(ierr, fnam) 359 ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam)352 ierr=NF90_INQ_VARID(ncid, trim(varname), varid); CALL ncerr(ierr, fnam) 360 353 ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam) 354 355 !--- Read unit for sea-ice variable only 356 IF (mode=='SIC') THEN 357 ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic) 358 IF(ierr/=NF90_NOERR) THEN 359 IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value' 360 unit_sic='X' 361 ELSE 362 IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic 363 END IF 364 END IF 361 365 362 366 !--- Longitude … … 365 369 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) 366 370 ierr=NF90_GET_VAR(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam) 367 WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep371 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep 368 372 369 373 !--- Latitude … … 372 376 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) 373 377 ierr=NF90_GET_VAR(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam) 374 WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep378 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep 375 379 376 380 !--- Time (variable is not needed - it is rebuilt - but calendar is) … … 385 389 CASE('SIC', 'SST'); cal_in='gregorian' 386 390 END SELECT 387 WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &391 IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & 388 392 // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.' 389 393 END IF 390 WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &394 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & 391 395 cal_in 392 396 397 393 398 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION -------------------- 394 399 !--- Determining input file number of days, depending on calendar … … 398 403 !--- If input records are not monthly, time sampling has to be constant ! 399 404 timeyear=mid_months(anneeref, cal_in, lmdep) 400 IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' & 401 // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, & 402 ' enregistrements.' 405 IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), & 406 ' ne comportent pas 12, mais ', lmdep, ' enregistrements.' 403 407 404 408 !--- GETTING THE FIELD AND INTERPOLATING IT ---------------------------------- … … 406 410 IF(extrp) ALLOCATE(work(imdep, jmdep)) 407 411 408 WRITE(lunout, *) 409 WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, & 410 ' CHAMPS.' 412 IF (prt_level>5) WRITE(lunout, *) 413 IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.' 411 414 ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) 412 415 DO l=1, lmdep … … 419 422 work) 420 423 421 IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN 422 IF(l==1) THEN 423 WRITE(lunout, *) & 424 '-------------------------------------------------------------------------' 425 WRITE(lunout, *) & 426 'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$' 427 WRITE(lunout, *) & 428 '-------------------------------------------------------------------------' 424 IF(ibar .AND. .NOT.oldice) THEN 425 IF(l==1 .AND. prt_level>5) THEN 426 WRITE(lunout, *) '-------------------------------------------------------------------------' 427 WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$' 428 WRITE(lunout, *) '-------------------------------------------------------------------------' 429 429 END IF 430 430 IF(mode=='RUG') champ=LOG(champ) … … 453 453 454 454 !--- TIME INTERPOLATION ------------------------------------------------------ 455 WRITE(lunout, *) 456 WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' 457 WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear 458 WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays 455 IF (prt_level>5) THEN 456 WRITE(lunout, *) 457 WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' 458 WRITE(lunout, *)' Vecteur temps en entree: ', timeyear 459 WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays 460 END IF 461 459 462 ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays)) 460 463 skip = .false. … … 471 474 END DO 472 475 if (n_extrap /= 0) then 473 print *,"get_2Dfield pchfe_95: n_extrap = ", n_extrap476 WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap 474 477 end if 475 478 champan(iip1, :, :)=champan(1, :, :) … … 479 482 DO j=1, jjp1 480 483 CALL minmax(iip1, champan(1, j, 10), chmin, chmax) 481 WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j484 IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j 482 485 END DO 483 486 484 487 !--- SPECIAL FILTER FOR SST: SST>271.38 -------------------------------------- 485 488 IF(mode=='SST') THEN 486 WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'489 IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38' 487 490 WHERE(champan<271.38) champan=271.38 488 491 END IF … … 490 493 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ------------------------------------- 491 494 IF(mode=='SIC') THEN 492 WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' 493 IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100. 495 IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' 496 497 IF (unit_sic=='1') THEN 498 ! Nothing to be done for sea-ice field is already in fraction of 1 499 ! This is the case for sea-ice in file cpl_atm_sic.nc 500 IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1' 501 ELSE 502 ! Convert sea ice from percentage to fraction of 1 503 IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' 504 champan(:, :, :)=champan(:, :, :)/100. 505 END IF 506 494 507 champan(iip1, :, :)=champan(1, :, :) 495 508 WHERE(champan>1.0) champan=1.0 496 509 WHERE(champan<0.0) champan=0.0 497 510 END IF 498 511 499 512 !--- DYNAMICAL TO PHYSICAL GRID ---------------------------------------------- -
trunk/LMDZ.COMMON/libf/dyn3d/vlsplt.F
r109 r270 1 ! 2 ! $Header$ 3 ! 4 c 1 c 2 c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z lguez $ 5 3 c 6 4 … … 131 129 IMPLICIT NONE 132 130 c 133 #include "dimensions.h" 134 #include "paramet.h" 135 #include "logic.h" 136 #include "comvert.h" 137 #include "comconst.h" 131 include "dimensions.h" 132 include "paramet.h" 133 include "logic.h" 134 include "comvert.h" 135 include "comconst.h" 136 include "iniprint.h" 138 137 c 139 138 c … … 353 352 354 353 IF(n0.gt.0) THEN 355 PRINT*,'Nombre de points pour lesquels on advect plus que le' 354 if (prt_level > 2) PRINT *, 355 $ 'Nombre de points pour lesquels on advect plus que le' 356 356 & ,'contenu de la maille : ',n0 357 357 -
trunk/LMDZ.COMMON/libf/dyn3d/write_paramLMDZ_dyn.h
r1 r270 51 51 . zx_tmp_2d,iip1*jjp1,ndex2d) 52 52 c 53 zx_tmp_2d(1:iip1,1:jjp1)=REAL( idissip)54 CALL histwrite(nid_ctesGCM, " idissip", itau_w,53 zx_tmp_2d(1:iip1,1:jjp1)=REAL(dissip_period) 54 CALL histwrite(nid_ctesGCM, "dissip_period", itau_w, 55 55 . zx_tmp_2d,iip1*jjp1,ndex2d) 56 56 c -
trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90
r269 r270 1 ! 2 ! $Id: advtrac_p.F 1446 2010-10-22 09:27:25Z emillour $ 3 ! 4 c 5 c 6 SUBROUTINE advtrac_p(pbaru,pbarv , 7 * p, masse,q,iapptrac,teta, 8 * flxw, 9 * pk ) 10 11 c Auteur : F. Hourdin 12 c 13 c Modif. P. Le Van (20/12/97) 14 c F. Codron (10/99) 15 c D. Le Croller (07/2001) 16 c M.A Filiberti (04/2002) 17 c 18 USE parallel 19 USE Write_Field_p 20 USE Bands 21 USE mod_hallo 22 USE Vampir 23 USE times 24 USE infotrac 25 USE control_mod 26 IMPLICIT NONE 27 c 28 #include "dimensions.h" 29 #include "paramet.h" 30 #include "comconst.h" 31 #include "comvert.h" 32 #include "comdissip.h" 33 #include "comgeom2.h" 34 #include "logic.h" 35 #include "temps.h" 36 #include "ener.h" 37 #include "description.h" 38 39 c------------------------------------------------------------------- 40 c Arguments 41 c------------------------------------------------------------------- 42 c Ajout PPM 43 c-------------------------------------------------------- 44 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) 45 c-------------------------------------------------------- 46 INTEGER iapptrac 47 REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 48 REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm) 49 REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm) 50 REAL pk(ip1jmp1,llm) 51 REAL :: flxw(ip1jmp1,llm) 52 53 c------------------------------------------------------------- 54 c Variables locales 55 c------------------------------------------------------------- 56 57 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) 58 REAL massem(ip1jmp1,llm),zdp(ip1jmp1) 59 REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 60 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu 61 INTEGER iadvtr 62 INTEGER ij,l,iq,iiq 63 REAL zdpmin, zdpmax 64 SAVE iadvtr, massem, pbaruc, pbarvc 65 DATA iadvtr/0/ 66 c$OMP THREADPRIVATE(iadvtr) 67 c---------------------------------------------------------- 68 c Rajouts pour PPM 69 c---------------------------------------------------------- 70 INTEGER indice,n 71 REAL dtbon ! Pas de temps adaptatif pour que CFL<1 72 REAL CFLmaxz,aaa,bbb ! CFL maximum 73 REAL psppm(iim,jjp1) ! pression au sol 74 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) 75 REAL qppm(iim*jjp1,llm,nqtot) 76 REAL fluxwppm(iim,jjp1,llm) 77 REAL apppm(llmp1), bpppm(llmp1) 78 LOGICAL dum,fill 79 DATA fill/.true./ 80 DATA dum/.true./ 81 REAL,SAVE :: finmasse(ip1jmp1,llm) 82 integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j 83 type(Request) :: Request_vanleer 84 REAL,SAVE :: p_tmp( ip1jmp1,llmp1 ) 85 REAL,SAVE :: teta_tmp(ip1jmp1,llm) 86 REAL,SAVE :: pk_tmp(ip1jmp1,llm) 87 88 ijb_u=ij_begin 89 ije_u=ij_end 90 91 ijb_v=ij_begin-iip1 92 ije_v=ij_end 93 if (pole_nord) ijb_v=ij_begin 94 if (pole_sud) ije_v=ij_end-iip1 95 96 IF(iadvtr.EQ.0) THEN 97 c CALL initial0(ijp1llm,pbaruc) 98 c CALL initial0(ijmllm,pbarvc) 99 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 100 DO l=1,llm 101 pbaruc(ijb_u:ije_u,l)=0. 102 pbarvc(ijb_v:ije_v,l)=0. 103 ENDDO 104 c$OMP END DO NOWAIT 105 ENDIF 106 107 c accumulation des flux de masse horizontaux 108 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 109 DO l=1,llm 110 DO ij = ijb_u,ije_u 111 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 112 ENDDO 113 DO ij = ijb_v,ije_v 114 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 115 ENDDO 116 ENDDO 117 c$OMP END DO NOWAIT 118 119 c selection de la masse instantannee des mailles avant le transport. 120 IF(iadvtr.EQ.0) THEN 121 122 c CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) 123 ijb=ij_begin 124 ije=ij_end 125 126 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 127 DO l=1,llm 128 massem(ijb:ije,l)=masse(ijb:ije,l) 129 ENDDO 130 c$OMP END DO NOWAIT 131 132 ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 133 c 134 ENDIF ! of IF(iadvtr.EQ.0) 135 136 iadvtr = iadvtr+1 137 138 c$OMP MASTER 139 iapptrac = iadvtr 140 c$OMP END MASTER 141 142 c Test pour savoir si on advecte a ce pas de temps 143 144 IF ( iadvtr.EQ.iapp_tracvl ) THEN 145 c$OMP MASTER 146 call suspend_timer(timer_caldyn) 147 c$OMP END MASTER 148 149 ijb=ij_begin 150 ije=ij_end 151 152 153 cc .. Modif P.Le Van ( 20/12/97 ) .... 154 cc 155 156 c traitement des flux de masse avant advection. 157 c 1. calcul de w 158 c 2. groupement des mailles pres du pole. 159 160 CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) 161 162 c$OMP BARRIER 163 164 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 165 DO l=1,llmp1 1 ! $Id: advtrac_p.F90 1549 2011-07-05 08:41:12Z lguez $ 2 3 SUBROUTINE advtrac_p(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk) 4 5 ! Auteur : F. Hourdin 6 ! 7 ! Modif. P. Le Van (20/12/97) 8 ! F. Codron (10/99) 9 ! D. Le Croller (07/2001) 10 ! M.A Filiberti (04/2002) 11 ! 12 USE parallel 13 USE Write_Field_p 14 USE Bands 15 USE mod_hallo 16 USE Vampir 17 USE times 18 USE infotrac 19 USE control_mod 20 IMPLICIT NONE 21 ! 22 include "dimensions.h" 23 include "paramet.h" 24 include "comconst.h" 25 include "comvert.h" 26 include "comdissip.h" 27 include "comgeom2.h" 28 include "logic.h" 29 include "temps.h" 30 include "ener.h" 31 include "description.h" 32 33 !------------------------------------------------------------------- 34 ! Arguments 35 !------------------------------------------------------------------- 36 ! Ajout PPM 37 !-------------------------------------------------------- 38 REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm) 39 !-------------------------------------------------------- 40 INTEGER iapptrac 41 REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) 42 REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm) 43 REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm) 44 REAL pk(ip1jmp1,llm) 45 REAL :: flxw(ip1jmp1,llm) 46 47 !------------------------------------------------------------- 48 ! Variables locales 49 !------------------------------------------------------------- 50 51 REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm) 52 REAL massem(ip1jmp1,llm),zdp(ip1jmp1) 53 REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm) 54 REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu 55 INTEGER iadvtr 56 INTEGER ij,l,iq,iiq 57 REAL zdpmin, zdpmax 58 SAVE iadvtr, massem, pbaruc, pbarvc 59 DATA iadvtr/0/ 60 !$OMP THREADPRIVATE(iadvtr) 61 !---------------------------------------------------------- 62 ! Rajouts pour PPM 63 !---------------------------------------------------------- 64 INTEGER indice,n 65 REAL dtbon ! Pas de temps adaptatif pour que CFL<1 66 REAL CFLmaxz,aaa,bbb ! CFL maximum 67 REAL psppm(iim,jjp1) ! pression au sol 68 REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm) 69 REAL qppm(iim*jjp1,llm,nqtot) 70 REAL fluxwppm(iim,jjp1,llm) 71 REAL apppm(llmp1), bpppm(llmp1) 72 LOGICAL dum,fill 73 DATA fill/.true./ 74 DATA dum/.true./ 75 REAL,SAVE :: finmasse(ip1jmp1,llm) 76 integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j 77 type(Request) :: Request_vanleer 78 REAL,SAVE :: p_tmp( ip1jmp1,llmp1 ) 79 REAL,SAVE :: teta_tmp(ip1jmp1,llm) 80 REAL,SAVE :: pk_tmp(ip1jmp1,llm) 81 82 ijb_u=ij_begin 83 ije_u=ij_end 84 85 ijb_v=ij_begin-iip1 86 ije_v=ij_end 87 if (pole_nord) ijb_v=ij_begin 88 if (pole_sud) ije_v=ij_end-iip1 89 90 IF(iadvtr.EQ.0) THEN 91 ! CALL initial0(ijp1llm,pbaruc) 92 ! CALL initial0(ijmllm,pbarvc) 93 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 94 DO l=1,llm 95 pbaruc(ijb_u:ije_u,l)=0. 96 pbarvc(ijb_v:ije_v,l)=0. 97 ENDDO 98 !$OMP END DO NOWAIT 99 ENDIF 100 101 ! accumulation des flux de masse horizontaux 102 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 103 DO l=1,llm 104 DO ij = ijb_u,ije_u 105 pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l) 106 ENDDO 107 DO ij = ijb_v,ije_v 108 pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l) 109 ENDDO 110 ENDDO 111 !$OMP END DO NOWAIT 112 113 ! selection de la masse instantannee des mailles avant le transport. 114 IF(iadvtr.EQ.0) THEN 115 116 ! CALL SCOPY(ip1jmp1*llm,masse,1,massem,1) 117 ijb=ij_begin 118 ije=ij_end 119 120 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 121 DO l=1,llm 122 massem(ijb:ije,l)=masse(ijb:ije,l) 123 ENDDO 124 !$OMP END DO NOWAIT 125 126 !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 ) 127 ! 128 ENDIF ! of IF(iadvtr.EQ.0) 129 130 iadvtr = iadvtr+1 131 132 !$OMP MASTER 133 iapptrac = iadvtr 134 !$OMP END MASTER 135 136 ! Test pour savoir si on advecte a ce pas de temps 137 138 IF ( iadvtr.EQ.iapp_tracvl ) THEN 139 !$OMP MASTER 140 call suspend_timer(timer_caldyn) 141 !$OMP END MASTER 142 143 ijb=ij_begin 144 ije=ij_end 145 146 147 !c .. Modif P.Le Van ( 20/12/97 ) .... 148 !c 149 150 ! traitement des flux de masse avant advection. 151 ! 1. calcul de w 152 ! 2. groupement des mailles pres du pole. 153 154 CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg ) 155 156 !$OMP BARRIER 157 158 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 159 DO l=1,llmp1 166 160 p_tmp(ijb:ije,l)=p(ijb:ije,l) 167 168 c$OMP END DO NOWAIT169 170 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)171 161 ENDDO 162 !$OMP END DO NOWAIT 163 164 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 165 DO l=1,llm 172 166 pk_tmp(ijb:ije,l)=pk(ijb:ije,l) 173 167 teta_tmp(ijb:ije,l)=teta(ijb:ije,l) 174 175 c$OMP END DO NOWAIT176 177 c$OMP MASTER178 179 c$OMP END MASTER180 181 call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,182 *jj_Nb_vanleer,0,0,Request_vanleer)183 call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,184 *jj_Nb_vanleer,1,0,Request_vanleer)185 call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,186 *jj_Nb_vanleer,0,0,Request_vanleer)187 call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,188 *jj_Nb_vanleer,0,0,Request_vanleer)189 call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,190 *jj_Nb_vanleer,1,1,Request_vanleer)191 call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,192 *jj_Nb_vanleer,1,1,Request_vanleer)193 call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,194 *jj_Nb_vanleer,1,1,Request_vanleer)195 196 call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, 197 *jj_nb_vanleer,0,0,Request_vanleer)198 199 200 201 c$OMP BARRIER202 203 204 205 c$OMP BARRIER206 c$OMP MASTER207 208 209 210 211 c$OMP END MASTER212 c$OMP BARRIER213 214 215 216 217 218 219 ctest sur l'eventuelle creation de valeurs negatives de la masse220 221 222 223 224 225 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)226 227 228 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)229 s - pbarvg(ij-iip1,l) + pbarvg(ij,l)230 s+ wg(ij,l+1) - wg(ij,l)231 232 233 cCALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )234 cym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?235 236 237 238 239 240 241 242 ENDDO243 244 245 cCALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )246 cym ---> eventuellement a revoir247 248 249 250 PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,251 s' MAX:', zdpmax252 253 254 255 c$OMP END DO NOWAIT256 257 c-------------------------------------------------------------------258 cAdvection proprement dite (Modification Le Croller (07/2001)259 c-------------------------------------------------------------------260 261 c----------------------------------------------------262 cCalcul des moyennes basées sur la masse263 c----------------------------------------------------264 265 cym ----> Normalement, inutile pour les schémas classiques266 cym ----> Revérifier lors de la parallélisation des autres schemas267 268 cym call massbar_p(massem,massebx,masseby)269 270 call vlspltgen_p( q,iadv, 2., massem, wg ,271 *pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )272 273 274 275 c-----------------------------------------------------------276 cAppel des sous programmes d'advection277 c-----------------------------------------------------------278 279 ccall clock(t_initial)168 ENDDO 169 !$OMP END DO NOWAIT 170 171 !$OMP MASTER 172 call VTb(VTHallo) 173 !$OMP END MASTER 174 175 call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm, & 176 jj_Nb_vanleer,0,0,Request_vanleer) 177 call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm, & 178 jj_Nb_vanleer,1,0,Request_vanleer) 179 call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm, & 180 jj_Nb_vanleer,0,0,Request_vanleer) 181 call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm, & 182 jj_Nb_vanleer,0,0,Request_vanleer) 183 call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm, & 184 jj_Nb_vanleer,1,1,Request_vanleer) 185 call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1, & 186 jj_Nb_vanleer,1,1,Request_vanleer) 187 call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm, & 188 jj_Nb_vanleer,1,1,Request_vanleer) 189 do j=1,nqtot 190 call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, & 191 jj_nb_vanleer,0,0,Request_vanleer) 192 enddo 193 194 call SendRequest(Request_vanleer) 195 !$OMP BARRIER 196 call WaitRequest(Request_vanleer) 197 198 199 !$OMP BARRIER 200 !$OMP MASTER 201 call SetDistrib(jj_nb_vanleer) 202 call VTe(VTHallo) 203 call VTb(VTadvection) 204 call start_timer(timer_vanleer) 205 !$OMP END MASTER 206 !$OMP BARRIER 207 208 ! ... Flux de masse diaganostiques traceurs 209 ijb=ij_begin 210 ije=ij_end 211 flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl) 212 213 ! test sur l'eventuelle creation de valeurs negatives de la masse 214 ijb=ij_begin 215 ije=ij_end 216 if (pole_nord) ijb=ij_begin+iip1 217 if (pole_sud) ije=ij_end-iip1 218 219 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 220 DO l=1,llm-1 221 DO ij = ijb+1,ije 222 zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) & 223 - pbarvg(ij-iip1,l) + pbarvg(ij,l) & 224 + wg(ij,l+1) - wg(ij,l) 225 ENDDO 226 227 ! CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 ) 228 ! ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ? 229 230 do ij=ijb,ije-iip1+1,iip1 231 zdp(ij)=zdp(ij+iip1-1) 232 enddo 233 234 DO ij = ijb,ije 235 zdp(ij)= zdp(ij)*dtvr/ massem(ij,l) 236 ENDDO 237 238 239 ! CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax ) 240 ! ym ---> eventuellement a revoir 241 CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax ) 242 243 IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN 244 PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, & 245 ' MAX:', zdpmax 246 ENDIF 247 248 ENDDO 249 !$OMP END DO NOWAIT 250 251 !------------------------------------------------------------------- 252 ! Advection proprement dite (Modification Le Croller (07/2001) 253 !------------------------------------------------------------------- 254 255 !---------------------------------------------------- 256 ! Calcul des moyennes basées sur la masse 257 !---------------------------------------------------- 258 259 !ym ----> Normalement, inutile pour les schémas classiques 260 !ym ----> Revérifier lors de la parallélisation des autres schemas 261 262 !ym call massbar_p(massem,massebx,masseby) 263 264 call vlspltgen_p( q,iadv, 2., massem, wg , & 265 pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp ) 266 267 268 GOTO 1234 269 !----------------------------------------------------------- 270 ! Appel des sous programmes d'advection 271 !----------------------------------------------------------- 272 do iq=1,nqtot 273 ! call clock(t_initial) 280 274 if(iadv(iq) == 0) cycle 281 c----------------------------------------------------------------282 cSchema de Van Leer I MUSCL283 c----------------------------------------------------------------275 ! ---------------------------------------------------------------- 276 ! Schema de Van Leer I MUSCL 277 ! ---------------------------------------------------------------- 284 278 if(iadv(iq).eq.10) THEN 285 286 287 288 c----------------------------------------------------------------289 cSchema "pseudo amont" + test sur humidite specifique290 Cpour la vapeur d'eau. F. Codron291 c----------------------------------------------------------------279 280 call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr) 281 282 ! ---------------------------------------------------------------- 283 ! Schema "pseudo amont" + test sur humidite specifique 284 ! pour la vapeur d'eau. F. Codron 285 ! ---------------------------------------------------------------- 292 286 else if(iadv(iq).eq.14) then 293 c 294 cymstop 'advtrac : appel à vlspltqs :schema non parallelise'295 CALL vlspltqs_p( q(1,1,1), 2., massem, wg , 296 *pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )297 c----------------------------------------------------------------298 cSchema de Frederic Hourdin299 c----------------------------------------------------------------287 ! 288 !ym stop 'advtrac : appel à vlspltqs :schema non parallelise' 289 CALL vlspltqs_p( q(1,1,1), 2., massem, wg , & 290 pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp ) 291 ! ---------------------------------------------------------------- 292 ! Schema de Frederic Hourdin 293 ! ---------------------------------------------------------------- 300 294 else if(iadv(iq).eq.12) then 301 stop 'advtrac : schema non parallelise'302 cPas de temps adaptatif295 stop 'advtrac : schema non parallelise' 296 ! Pas de temps adaptatif 303 297 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 304 298 if (n.GT.1) then 305 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',306 sdtvr,'n=',n299 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 300 dtvr,'n=',n 307 301 endif 308 302 do indice=1,n 309 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)303 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1) 310 304 end do 311 305 else if(iadv(iq).eq.13) then 312 stop 'advtrac : schema non parallelise'313 cPas de temps adaptatif306 stop 'advtrac : schema non parallelise' 307 ! Pas de temps adaptatif 314 308 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 315 309 if (n.GT.1) then 316 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',317 sdtvr,'n=',n310 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 311 dtvr,'n=',n 318 312 endif 319 do indice=1,n320 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)321 end do322 c----------------------------------------------------------------323 cSchema de pente SLOPES324 c----------------------------------------------------------------313 do indice=1,n 314 call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2) 315 end do 316 ! ---------------------------------------------------------------- 317 ! Schema de pente SLOPES 318 ! ---------------------------------------------------------------- 325 319 else if (iadv(iq).eq.20) then 326 stop 'advtrac : schema non parallelise'327 328 329 330 c----------------------------------------------------------------331 cSchema de Prather332 c----------------------------------------------------------------320 stop 'advtrac : schema non parallelise' 321 322 call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0) 323 324 ! ---------------------------------------------------------------- 325 ! Schema de Prather 326 ! ---------------------------------------------------------------- 333 327 else if (iadv(iq).eq.30) then 334 stop 'advtrac : schema non parallelise'335 cPas de temps adaptatif328 stop 'advtrac : schema non parallelise' 329 ! Pas de temps adaptatif 336 330 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 337 331 if (n.GT.1) then 338 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',339 sdtvr,'n=',n332 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 333 dtvr,'n=',n 340 334 endif 341 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, 342 sn,dtbon)343 c----------------------------------------------------------------344 cSchemas PPM Lin et Rood345 c----------------------------------------------------------------346 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.347 siadv(iq).LE.18)) then335 call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, & 336 n,dtbon) 337 ! ---------------------------------------------------------------- 338 ! Schemas PPM Lin et Rood 339 ! ---------------------------------------------------------------- 340 else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. & 341 iadv(iq).LE.18)) then 348 342 349 343 stop 'advtrac : schema non parallelise' 350 344 351 c Test sur le flux horizontal 352 c Pas de temps adaptatif 353 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 354 if (n.GT.1) then 355 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', 356 s dtvr,'n=',n 357 endif 358 c Test sur le flux vertical 359 CFLmaxz=0. 360 do l=2,llm 361 do ij=iip2,ip1jm 362 aaa=wg(ij,l)*dtvr/massem(ij,l) 363 CFLmaxz=max(CFLmaxz,aaa) 364 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 365 CFLmaxz=max(CFLmaxz,bbb) 345 ! Test sur le flux horizontal 346 ! Pas de temps adaptatif 347 call adaptdt(iadv(iq),dtbon,n,pbarug,massem) 348 if (n.GT.1) then 349 write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', & 350 dtvr,'n=',n 351 endif 352 ! Test sur le flux vertical 353 CFLmaxz=0. 354 do l=2,llm 355 do ij=iip2,ip1jm 356 aaa=wg(ij,l)*dtvr/massem(ij,l) 357 CFLmaxz=max(CFLmaxz,aaa) 358 bbb=-wg(ij,l)*dtvr/massem(ij,l-1) 359 CFLmaxz=max(CFLmaxz,bbb) 360 enddo 366 361 enddo 367 enddo 368 if (CFLmaxz.GE.1) then 369 write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 370 endif 371 372 c----------------------------------------------------------- 373 c Ss-prg interface LMDZ.4->PPM3d 374 c----------------------------------------------------------- 375 376 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, 377 s apppm,bpppm,massebx,masseby,pbarug,pbarvg, 378 s unatppm,vnatppm,psppm) 379 380 do indice=1,n 381 c--------------------------------------------------------------------- 382 c VL (version PPM) horiz. et PPM vert. 383 c--------------------------------------------------------------------- 384 if (iadv(iq).eq.11) then 385 c Ss-prg PPM3d de Lin 386 call ppm3d(1,qppm(1,1,iq), 387 s psppm,psppm, 388 s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, 389 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 390 s fill,dum,220.) 391 392 c---------------------------------------------------------------------- 393 c Monotonic PPM 394 c---------------------------------------------------------------------- 395 else if (iadv(iq).eq.16) then 396 c Ss-prg PPM3d de Lin 397 call ppm3d(1,qppm(1,1,iq), 398 s psppm,psppm, 399 s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, 400 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 401 s fill,dum,220.) 402 c--------------------------------------------------------------------- 403 404 c--------------------------------------------------------------------- 405 c Semi Monotonic PPM 406 c--------------------------------------------------------------------- 407 else if (iadv(iq).eq.17) then 408 c Ss-prg PPM3d de Lin 409 call ppm3d(1,qppm(1,1,iq), 410 s psppm,psppm, 411 s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, 412 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 413 s fill,dum,220.) 414 c--------------------------------------------------------------------- 415 416 c--------------------------------------------------------------------- 417 c Positive Definite PPM 418 c--------------------------------------------------------------------- 419 else if (iadv(iq).eq.18) then 420 c Ss-prg PPM3d de Lin 421 call ppm3d(1,qppm(1,1,iq), 422 s psppm,psppm, 423 s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, 424 s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, 425 s fill,dum,220.) 426 c--------------------------------------------------------------------- 427 endif 428 enddo 429 c----------------------------------------------------------------- 430 c Ss-prg interface PPM3d-LMDZ.4 431 c----------------------------------------------------------------- 432 call interpost(q(1,1,iq),qppm(1,1,iq)) 433 endif 434 c---------------------------------------------------------------------- 435 436 c----------------------------------------------------------------- 437 c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 438 c et Nord j=1 439 c----------------------------------------------------------------- 440 441 c call traceurpole(q(1,1,iq),massem) 442 443 c calcul du temps cpu pour un schema donne 444 445 c call clock(t_final) 446 cym tps_cpu=t_final-t_initial 447 cym cpuadv(iq)=cpuadv(iq)+tps_cpu 448 449 end DO 450 451 1234 CONTINUE 452 c$OMP BARRIER 453 454 if (planet_type=="earth") then 362 if (CFLmaxz.GE.1) then 363 write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz 364 endif 365 366 !----------------------------------------------------------- 367 ! Ss-prg interface LMDZ.4->PPM3d 368 !----------------------------------------------------------- 369 370 call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, & 371 apppm,bpppm,massebx,masseby,pbarug,pbarvg, & 372 unatppm,vnatppm,psppm) 373 374 do indice=1,n 375 !---------------------------------------------------------------- 376 ! VL (version PPM) horiz. et PPM vert. 377 !---------------------------------------------------------------- 378 if (iadv(iq).eq.11) then 379 ! Ss-prg PPM3d de Lin 380 call ppm3d(1,qppm(1,1,iq), & 381 psppm,psppm, & 382 unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, & 383 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 384 fill,dum,220.) 385 386 !------------------------------------------------------------- 387 ! Monotonic PPM 388 !------------------------------------------------------------- 389 else if (iadv(iq).eq.16) then 390 ! Ss-prg PPM3d de Lin 391 call ppm3d(1,qppm(1,1,iq), & 392 psppm,psppm, & 393 unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, & 394 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 395 fill,dum,220.) 396 !------------------------------------------------------------- 397 398 !------------------------------------------------------------- 399 ! Semi Monotonic PPM 400 !------------------------------------------------------------- 401 else if (iadv(iq).eq.17) then 402 ! Ss-prg PPM3d de Lin 403 call ppm3d(1,qppm(1,1,iq), & 404 psppm,psppm, & 405 unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, & 406 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 407 fill,dum,220.) 408 !------------------------------------------------------------- 409 410 !------------------------------------------------------------- 411 ! Positive Definite PPM 412 !------------------------------------------------------------- 413 else if (iadv(iq).eq.18) then 414 ! Ss-prg PPM3d de Lin 415 call ppm3d(1,qppm(1,1,iq), & 416 psppm,psppm, & 417 unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, & 418 iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, & 419 fill,dum,220.) 420 !------------------------------------------------------------- 421 endif 422 enddo 423 !----------------------------------------------------------------- 424 ! Ss-prg interface PPM3d-LMDZ.4 425 !----------------------------------------------------------------- 426 call interpost(q(1,1,iq),qppm(1,1,iq)) 427 endif 428 !---------------------------------------------------------------------- 429 430 !----------------------------------------------------------------- 431 ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 432 ! et Nord j=1 433 !----------------------------------------------------------------- 434 435 ! call traceurpole(q(1,1,iq),massem) 436 437 ! calcul du temps cpu pour un schema donne 438 439 ! call clock(t_final) 440 !ym tps_cpu=t_final-t_initial 441 !ym cpuadv(iq)=cpuadv(iq)+tps_cpu 442 443 end DO 444 445 1234 CONTINUE 446 !$OMP BARRIER 447 448 if (planet_type=="earth") then 455 449 456 450 ijb=ij_begin 457 451 ije=ij_end 458 452 459 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)453 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 460 454 DO l = 1, llm 461 DO ij = ijb, ije462 finmasse(ij,l) = p(ij,l) - p(ij,l+1)463 ENDDO455 DO ij = ijb, ije 456 finmasse(ij,l) = p(ij,l) - p(ij,l+1) 457 ENDDO 464 458 ENDDO 465 c$OMP END DO459 !$OMP END DO 466 460 467 461 CALL qminimum_p( q, 2, finmasse ) 468 462 469 c------------------------------------------------------------------470 con reinitialise a zero les flux de masse cumules471 c---------------------------------------------------472 ciadvtr=0473 474 c$OMP MASTER463 !------------------------------------------------------------------ 464 ! on reinitialise a zero les flux de masse cumules 465 !--------------------------------------------------- 466 ! iadvtr=0 467 468 !$OMP MASTER 475 469 call VTe(VTadvection) 476 470 call stop_timer(timer_vanleer) 477 471 call VTb(VThallo) 478 c$OMP END MASTER472 !$OMP END MASTER 479 473 480 474 do j=1,nqtot 481 call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,482 *jj_nb_caldyn,0,0,Request_vanleer)475 call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, & 476 jj_nb_caldyn,0,0,Request_vanleer) 483 477 enddo 484 478 485 call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, 486 *jj_nb_caldyn,0,0,Request_vanleer)479 call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, & 480 jj_nb_caldyn,0,0,Request_vanleer) 487 481 488 482 call SendRequest(Request_vanleer) 489 c$OMP BARRIER483 !$OMP BARRIER 490 484 call WaitRequest(Request_vanleer) 491 485 492 c$OMP BARRIER493 c$OMP MASTER486 !$OMP BARRIER 487 !$OMP MASTER 494 488 call SetDistrib(jj_nb_caldyn) 495 489 call VTe(VThallo) 496 490 call resume_timer(timer_caldyn) 497 c$OMP END MASTER 498 c$OMP BARRIER 499 iadvtr=0 500 endif ! of if (planet_type=="earth") 501 ENDIF ! if iadvtr.EQ.iapp_tracvl 502 503 RETURN 504 END 505 491 !$OMP END MASTER 492 !$OMP BARRIER 493 iadvtr=0 494 endif ! of if (planet_type=="earth") 495 ENDIF ! if iadvtr.EQ.iapp_tracvl 496 497 END SUBROUTINE advtrac_p -
trunk/LMDZ.COMMON/libf/dyn3dpar/bilan_dyn_p.F
r1 r270 389 389 Q_cum(:,jjb:jje,:,:)=0. 390 390 flux_uQ_cum(:,jjb:jje,:,:)=0. 391 flux_v_cum(:,jjb:jje,:)=0.392 391 if (pole_sud) jje=jj_end-1 393 392 flux_v_cum(:,jjb:jje,:)=0. -
trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90
r101 r270 1 1 ! 2 ! $Id: ce0l.F90 1 482 2011-02-09 15:03:09Z jghattas $2 ! $Id: ce0l.F90 1511 2011-04-28 15:21:47Z jghattas $ 3 3 ! 4 4 !------------------------------------------------------------------------------- … … 19 19 USE dimphy 20 20 USE comgeomphy 21 USE mod_phys_lmdz_para22 USE mod_const_mpi23 21 USE infotrac 24 USE parallel, ONLY: finalize_parallel25 22 26 23 #ifdef CPP_IOIPSL … … 30 27 #endif 31 28 IMPLICIT NONE 32 #include "iniprint.h"33 29 #ifndef CPP_EARTH 34 30 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics' … … 40 36 #include "paramet.h" 41 37 #include "indicesol.h" 38 #include "iniprint.h" 42 39 #include "temps.h" 43 40 #include "logic.h" 41 INTEGER, PARAMETER :: longcles=20 42 REAL, DIMENSION(longcles) :: clesphy0 44 43 REAL, DIMENSION(iip1,jjp1) :: masque 45 44 CHARACTER(LEN=15) :: calnd 45 REAL, DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol 46 46 !------------------------------------------------------------------------------- 47 CALL conf_gcm( 99, .TRUE. ) 48 49 CALL init_mpi 47 CALL conf_gcm( 99, .TRUE. , clesphy0 ) 50 48 51 49 CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/)) 52 50 WRITE(lunout,*)'---> klon=',klon 53 IF (mpi_size>1 .OR. omp_size>1) THEN54 CALL abort_gcm('ce0l','In parallel mode, &55 & ce0l must be called only &56 & for 1 process and 1 task',1)57 ENDIF58 59 51 CALL InitComgeomphy 60 52 … … 89 81 WRITE(lunout,'(//)') 90 82 WRITE(lunout,*) ' interbar = ',interbar 91 CALL etat0_netcdf(interbar,masque, ok_etat0)83 CALL etat0_netcdf(interbar,masque,phis,ok_etat0) 92 84 93 85 IF(ok_limit) THEN … … 100 92 END IF 101 93 102 !$OMP MASTER 103 CALL finalize_parallel 104 !$OMP END MASTER 105 94 IF (grilles_gcm_netcdf) THEN 95 WRITE(lunout,'(//)') 96 WRITE(lunout,*) ' *************************** ' 97 WRITE(lunout,*) ' *** grilles_gcm_netcdf *** ' 98 WRITE(lunout,*) ' *************************** ' 99 WRITE(lunout,'(//)') 100 CALL grilles_gcm_netcdf_sub(masque,phis) 101 END IF 106 102 #endif 107 103 ! of #ifndef CPP_EARTH #else -
trunk/LMDZ.COMMON/libf/dyn3dpar/comdissipn.h
r1 r270 2 2 ! $Header$ 3 3 ! 4 c----------------------------------------------------------------------- 5 c INCLUDE comdissipn.h 4 ! Attention : ce fichier include est compatible format fixe/format libre 5 ! veillez à n'utiliser que des ! pour les commentaires 6 ! et à bien positionner les & des lignes de continuation 7 ! (les placer en colonne 6 et en colonne 73) 8 !----------------------------------------------------------------------- 9 ! INCLUDE comdissipn.h 6 10 7 11 REAL tetaudiv, tetaurot, tetah, cdivu, crot, cdivh 8 c 9 COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm) , 10 1cdivu, crot, cdivh12 ! 13 COMMON/comdissipn/ tetaudiv(llm),tetaurot(llm),tetah(llm) , & 14 & cdivu, crot, cdivh 11 15 12 c 13 cLes parametres de ce common proviennent des calculs effectues dans14 cInidissip .15 c 16 c-----------------------------------------------------------------------16 ! 17 ! Les parametres de ce common proviennent des calculs effectues dans 18 ! Inidissip . 19 ! 20 !----------------------------------------------------------------------- -
trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
r119 r270 1 1 ! 2 ! $Id: conf_gcm.F 14 38 2010-10-08 10:19:34Z jghattas $2 ! $Id: conf_gcm.F 1418 2010-07-19 15:11:24Z jghattas $ 3 3 ! 4 4 c … … 6 6 SUBROUTINE conf_gcm( tapedef, etatinit ) 7 7 c 8 USE control_mod 8 9 #ifdef CPP_IOIPSL 9 10 use IOIPSL … … 12 13 use ioipsl_getincom 13 14 #endif 14 use misc_mod15 use mod_filtre_fft, ONLY : use_filtre_fft16 use mod_hallo, ONLY : use_mpi_alloc17 use parallel, ONLY : omp_chunk18 USE control_mod19 15 IMPLICIT NONE 20 16 c----------------------------------------------------------------------- … … 37 33 #include "serre.h" 38 34 #include "comdissnew.h" 39 !#include "clesphys.h"40 #include "iniprint.h"41 35 #include "temps.h" 42 36 #include "comconst.h" 43 37 44 38 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique 39 ! #include "clesphys.h" 40 #include "iniprint.h" 45 41 c 46 42 c … … 53 49 LOGICAL fxyhypbb, ysinuss 54 50 INTEGER i 55 51 LOGICAL use_filtre_fft 56 52 c 57 53 c ------------------------------------------------------------------- … … 82 78 c initialisations: 83 79 c ---------------- 84 adjust=.false. 85 call getin('adjust',adjust) 86 87 itaumax=0 88 call getin('itaumax',itaumax); 89 if (itaumax<=0) itaumax=HUGE(itaumax) 90 80 91 81 !Config Key = lunout 92 82 !Config Desc = unite de fichier pour les impressions … … 190 180 191 181 !Config Key = nsplit_phys 192 !Config Desc = nombre d'iteration de la physique 193 !Config Def = 240 194 !Config Help = nombre d'itration de la physique 195 ! 182 !Config Desc = nombre de subdivisions par pas physique 183 !Config Def = 1 184 !Config Help = nombre de subdivisions par pas physique 196 185 nsplit_phys = 1 197 186 CALL getin('nsplit_phys',nsplit_phys) … … 241 230 CALL getin('output_grads_dyn',output_grads_dyn) 242 231 243 !Config Key = idissip232 !Config Key = dissip_period 244 233 !Config Desc = periode de la dissipation 245 !Config Def = 10234 !Config Def = 0 246 235 !Config Help = periode de la dissipation 247 !Config (en pas) ... a completer ! 248 idissip = 10 249 CALL getin('idissip',idissip) 236 !Config dissip_period=0 => la valeur sera calcule dans inidissip 237 !Config dissip_period>0 => on prend cette valeur 238 dissip_period = 0 239 CALL getin('dissip_period',dissip_period) 250 240 251 241 ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... … … 340 330 CALL getin('tau_top_bound',tau_top_bound) 341 331 342 !343 332 !Config Key = coefdis 344 333 !Config Desc = coefficient pour gamdissip … … 407 396 ip_ebil_dyn = 0 408 397 CALL getin('ip_ebil_dyn',ip_ebil_dyn) 409 !410 411 398 412 399 ccc .... P. Le Van , ajout le 7/03/95 .pour le zoom ... … … 613 600 offline = .FALSE. 614 601 CALL getin('offline',offline) 615 IF (offline .AND. adjust) THEN 616 WRITE(lunout,*) 617 & 'WARNING : option offline does not work with adjust=y :' 618 WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ', 619 & 'and fluxstokev.nc will not be created' 620 WRITE(lunout,*) 621 & 'only the file phystoke.nc will still be created ' 622 END IF 623 602 624 603 !Config Key = config_inca 625 604 !Config Desc = Choix de configuration de INCA … … 655 634 ok_dyn_ave = .FALSE. 656 635 CALL getin('ok_dyn_ave',ok_dyn_ave) 636 657 637 658 638 write(lunout,*)' #########################################' … … 669 649 write(lunout,*)' day_step = ', day_step 670 650 write(lunout,*)' iperiod = ', iperiod 671 write(lunout,*)' nsplit_phys = ', nsplit_phys672 651 write(lunout,*)' iconser = ', iconser 673 652 write(lunout,*)' iecri = ', iecri 674 653 write(lunout,*)' periodav = ', periodav 675 654 write(lunout,*)' output_grads_dyn = ', output_grads_dyn 676 write(lunout,*)' idissip = ', idissip655 write(lunout,*)' dissip_period = ', dissip_period 677 656 write(lunout,*)' lstardis = ', lstardis 678 657 write(lunout,*)' nitergdiv = ', nitergdiv … … 816 795 offline = .FALSE. 817 796 CALL getin('offline',offline) 818 IF (offline .AND. adjust) THEN819 WRITE(lunout,*)820 & 'WARNING : option offline does not work with adjust=y :'821 WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ',822 & 'and fluxstokev.nc will not be created'823 WRITE(lunout,*)824 & 'only the file phystoke.nc will still be created '825 END IF826 797 827 798 !Config Key = config_inca … … 836 807 837 808 !Config Key = ok_dynzon 838 !Config Desc = calcul et sortie des transports809 !Config Desc = sortie des transports zonaux dans la dynamique 839 810 !Config Def = n 840 !Config Help = Permet de mettre en route le calcul des transports811 !Config Help = 841 812 !Config 842 ok_dynzon = .FALSE.843 CALL getin('ok_dynzon',ok_dynzon)813 ok_dynzon = .FALSE. 814 CALL getin('ok_dynzon',ok_dynzon) 844 815 845 816 !Config Key = ok_dyn_ins … … 863 834 !Config Def = false 864 835 !Config Help = permet d'activer l'utilisation des FFT pour effectuer 865 !Config le filtrage aux poles. 836 !Config le filtrage aux poles. 837 ! Le filtre fft n'est pas implemente dans dyn3d 866 838 use_filtre_fft=.FALSE. 867 839 CALL getin('use_filtre_fft',use_filtre_fft) 868 840 869 IF (use_filtre_fft .AND. grossismx /= 1.0) THEN 870 write(lunout,*)'WARNING !!! ' 871 write(lunout,*)"Le zoom en longitude est incompatible", 872 & " avec l'utilisation du filtre FFT ", 873 & "---> filtre FFT désactivé " 874 use_filtre_fft=.FALSE. 841 IF (use_filtre_fft) THEN 842 write(lunout,*)'STOP !!!' 843 write(lunout,*)'use_filtre_fft n est pas implemente dans dyn3d' 844 STOP 875 845 ENDIF 876 846 877 878 879 !Config Key = use_mpi_alloc880 !Config Desc = Utilise un buffer MPI en m�moire globale881 !Config Def = false882 !Config Help = permet d'activer l'utilisation d'un buffer MPI883 !Config en m�moire globale a l'aide de la fonction MPI_ALLOC.884 !Config Cela peut am�liorer la bande passante des transferts MPI885 !Config d'un facteur 2886 use_mpi_alloc=.FALSE.887 CALL getin('use_mpi_alloc',use_mpi_alloc)888 889 !Config Key = omp_chunk890 !Config Desc = taille des blocs openmp891 !Config Def = 1892 !Config Help = defini la taille des packets d'it�ration openmp893 !Config distribu�e � chaque t�che lors de l'entr�e dans une894 !Config boucle parall�lis�e895 896 omp_chunk=1897 CALL getin('omp_chunk',omp_chunk)898 899 847 !Config key = ok_strato 900 848 !Config Desc = activation de la version strato … … 902 850 !Config Help = active la version stratosphérique de LMDZ de F. Lott 903 851 904 ok_strato=. FALSE.852 ok_strato=.TRUE. 905 853 CALL getin('ok_strato',ok_strato) 906 854 … … 928 876 ok_etat0 = .TRUE. 929 877 CALL getin('ok_etat0',ok_etat0) 878 879 !Config Key = grilles_gcm_netcdf 880 !Config Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit 881 !Config Def = n 882 grilles_gcm_netcdf = .FALSE. 883 CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf) 930 884 931 885 write(lunout,*)' #########################################' … … 943 897 write(lunout,*)' periodav = ', periodav 944 898 write(lunout,*)' output_grads_dyn = ', output_grads_dyn 945 write(lunout,*)' idissip = ', idissip899 write(lunout,*)' dissip_period = ', dissip_period 946 900 write(lunout,*)' lstardis = ', lstardis 947 901 write(lunout,*)' nitergdiv = ', nitergdiv … … 968 922 write(lunout,*)' offline = ', offline 969 923 write(lunout,*)' config_inca = ', config_inca 970 write(lunout,*)' ok_dynzon = ', ok_dynzon 924 write(lunout,*)' ok_dynzon = ', ok_dynzon 971 925 write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins 972 926 write(lunout,*)' ok_dyn_ave = ', ok_dyn_ave 973 write(lunout,*)' use_filtre_fft = ', use_filtre_fft974 write(lunout,*)' use_mpi_alloc = ', use_mpi_alloc975 write(lunout,*)' omp_chunk = ', omp_chunk976 927 write(lunout,*)' ok_strato = ', ok_strato 977 928 write(lunout,*)' ok_gradsfile = ', ok_gradsfile 978 929 write(lunout,*)' ok_limit = ', ok_limit 979 930 write(lunout,*)' ok_etat0 = ', ok_etat0 931 write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf 980 932 c 981 933 RETURN -
trunk/LMDZ.COMMON/libf/dyn3dpar/control_mod.F90
r97 r270 12 12 REAL :: periodav 13 13 INTEGER :: nday,day_step,iperiod,iapp_tracvl,nsplit_phys 14 INTEGER :: iconser,iecri, idissip,iphysiq,iecrimoy14 INTEGER :: iconser,iecri,dissip_period,iphysiq,iecrimoy 15 15 INTEGER :: dayref,anneeref, raz_date, ip_ebil_dyn 16 16 LOGICAL :: offline -
trunk/LMDZ.COMMON/libf/dyn3dpar/defrun.F
r1 r270 6 6 SUBROUTINE defrun( tapedef, etatinit, clesphy0 ) 7 7 c 8 ! ========================== ATTENTION ============================= 9 ! COMMENTAIRE SL : 10 ! NE SERT PLUS APPAREMMENT 11 ! DONC PAS MIS A JOUR POUR L'UTILISATION AVEC LES PLANETES 12 ! ================================================================== 13 8 14 USE control_mod 15 9 16 IMPLICIT NONE 10 17 c----------------------------------------------------------------------- … … 132 139 133 140 READ (tapedef,9001) ch1,ch4 134 READ (tapedef,*) idissip135 WRITE(tapeout,9001) ch1,' idissip'136 WRITE(tapeout,*) idissip141 READ (tapedef,*) dissip_period 142 WRITE(tapeout,9001) ch1,'dissip_period' 143 WRITE(tapeout,*) dissip_period 137 144 138 145 ccc .... P. Le Van , modif le 29/04/97 .pour la dissipation ... … … 196 203 197 204 198 ccc .... P.Le Van, ajout le 03/01/96 pour l'ecriture phys ...199 c200 205 READ (tapedef,9001) ch1,ch4 201 206 READ (tapedef,*) cycle_diurne -
trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F
r127 r270 125 125 endif 126 126 !$OMP END MASTER 127 127 !$OMP BARRIER 128 128 jjb=jj_begin 129 129 jje=jj_end … … 171 171 endif 172 172 c$OMP END MASTER 173 c$OMP BARRIER 173 174 c 174 175 c -
trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F
r127 r270 121 121 endif 122 122 !$OMP END MASTER 123 123 !$OMP BARRIER 124 124 jjb=jj_begin 125 125 jje=jj_end … … 169 169 endif 170 170 c$OMP END MASTER 171 c$OMP BARRIER 171 172 c 172 173 c -
trunk/LMDZ.COMMON/libf/dyn3dpar/iniacademic.F90
r127 r270 1 1 ! 2 ! $Id: iniacademic.F90 152 0 2011-05-23 11:37:09Z emillour$2 ! $Id: iniacademic.F90 1529 2011-05-26 15:17:33Z fairhead $ 3 3 ! 4 4 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) … … 115 115 endif 116 116 117 academic_case: if (iflag_phys == 2) then117 academic_case: if (iflag_phys >= 2) then 118 118 ! initializations 119 119 … … 208 208 IF (.NOT. read_start) THEN 209 209 ! surface pressure 210 ps(:)=preff 210 if (iflag_phys>2) then 211 ps(:)=108080. ! Earth aqua/terra planets 212 else 213 ps(:)=preff 214 endif 211 215 ! ground geopotential 212 216 phis(:)=0. -
trunk/LMDZ.COMMON/libf/dyn3dpar/iniconst.F
r127 r270 47 47 c----------------------------------------------------------------------- 48 48 49 dtdiss = idissip * dtvr50 49 dtphys = iphysiq * dtvr 51 50 unsim = 1./iim -
trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90
r269 r270 1 1 ! 2 ! $Id: inidissip.F 1403 2010-07-01 09:02:53Z fairhead$2 ! $Id: inidissip.F90 1502 2011-03-21 16:07:54Z jghattas $ 3 3 ! 4 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , 5 * tetagdiv,tetagrot,tetatemp ) 6 c======================================================================= 7 c initialisation de la dissipation horizontale 8 c======================================================================= 9 c----------------------------------------------------------------------- 10 c declarations: 11 c ------------- 12 13 USE control_mod 14 15 IMPLICIT NONE 16 #include "dimensions.h" 17 #include "paramet.h" 18 #include "comdissipn.h" 19 #include "comconst.h" 20 #include "comvert.h" 21 #include "logic.h" 22 23 LOGICAL lstardis 24 INTEGER nitergdiv,nitergrot,niterh 25 REAL tetagdiv,tetagrot,tetatemp 26 REAL fact,zvert(llm),zz 27 REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm) 28 REAL ullm,vllm,umin,vmin,zhmin,zhmax 29 REAL zllm,z1llm 30 31 INTEGER l,ij,idum,ii 32 REAL tetamin 33 REAL Pup 34 35 REAL ran1 36 37 38 c----------------------------------------------------------------------- 39 c 40 c calcul des valeurs propres des operateurs par methode iterrative: 41 c ----------------------------------------------------------------- 42 43 crot = -1. 44 cdivu = -1. 45 cdivh = -1. 46 47 c calcul de la valeur propre de divgrad: 48 c -------------------------------------- 49 idum = 0 50 DO l = 1, llm 51 DO ij = 1, ip1jmp1 4 SUBROUTINE inidissip ( lstardis,nitergdiv,nitergrot,niterh , & 5 tetagdiv,tetagrot,tetatemp ) 6 !======================================================================= 7 ! initialisation de la dissipation horizontale 8 !======================================================================= 9 !----------------------------------------------------------------------- 10 ! declarations: 11 ! ------------- 12 13 USE control_mod, only : dissip_period,iperiod 14 15 IMPLICIT NONE 16 include "dimensions.h" 17 include "paramet.h" 18 include "comdissipn.h" 19 include "comconst.h" 20 include "comvert.h" 21 include "logic.h" 22 include "iniprint.h" 23 24 LOGICAL,INTENT(in) :: lstardis 25 INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh 26 REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp 27 28 ! Local variables: 29 REAL fact,zvert(llm),zz 30 REAL zh(ip1jmp1),zu(ip1jmp1),zv(ip1jm),deltap(ip1jmp1,llm) 31 REAL ullm,vllm,umin,vmin,zhmin,zhmax 32 REAL zllm,z1llm 33 34 INTEGER l,ij,idum,ii 35 REAL tetamin 36 REAL Pup 37 character (len=80) :: abort_message 38 39 REAL ran1 40 41 42 !----------------------------------------------------------------------- 43 ! 44 ! calcul des valeurs propres des operateurs par methode iterrative: 45 ! ----------------------------------------------------------------- 46 47 crot = -1. 48 cdivu = -1. 49 cdivh = -1. 50 51 ! calcul de la valeur propre de divgrad: 52 ! -------------------------------------- 53 idum = 0 54 DO l = 1, llm 55 DO ij = 1, ip1jmp1 52 56 deltap(ij,l) = 1. 53 ENDDO 54 ENDDO 55 56 idum = -1 57 zh(1) = RAN1(idum)-.5 58 idum = 0 59 DO ij = 2, ip1jmp1 60 zh(ij) = RAN1(idum) -.5 61 ENDDO 62 63 CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) 64 65 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 66 67 IF ( zhmin .GE. zhmax ) THEN 68 PRINT*,' Inidissip zh min max ',zhmin,zhmax 69 STOP'probleme generateur alleatoire dans inidissip' 70 ENDIF 71 72 zllm = ABS( zhmax ) 73 DO l = 1,50 74 IF(lstardis) THEN 75 CALL divgrad2(1,zh,deltap,niterh,zh) 76 ELSE 77 CALL divgrad (1,zh,niterh,zh) 78 ENDIF 79 80 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 81 82 zllm = ABS( zhmax ) 83 z1llm = 1./zllm 84 DO ij = 1,ip1jmp1 85 zh(ij) = zh(ij)* z1llm 86 ENDDO 87 ENDDO 88 89 IF(lstardis) THEN 90 cdivh = 1./ zllm 91 ELSE 92 cdivh = zllm ** ( -1./niterh ) 93 ENDIF 94 95 c calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 96 c ----------------------------------------------------------------- 97 print*,'calcul des valeurs propres' 98 99 DO 20 ii = 1, 2 100 c 101 DO ij = 1, ip1jmp1 102 zu(ij) = RAN1(idum) -.5 103 ENDDO 104 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 105 DO ij = 1, ip1jm 106 zv(ij) = RAN1(idum) -.5 107 ENDDO 108 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 109 110 CALL minmax(iip1*jjp1,zu,umin,ullm ) 111 CALL minmax(iip1*jjm, zv,vmin,vllm ) 112 113 ullm = ABS ( ullm ) 114 vllm = ABS ( vllm ) 115 116 DO 5 l = 1, 50 117 IF(ii.EQ.1) THEN 118 ccccc CALL covcont( 1,zu,zv,zu,zv ) 119 IF(lstardis) THEN 120 CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv ) 121 ELSE 122 CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv ) 123 ENDIF 124 ELSE 125 IF(lstardis) THEN 126 CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv ) 127 ELSE 128 CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv ) 129 ENDIF 130 ENDIF 131 132 CALL minmax(iip1*jjp1,zu,umin,ullm ) 133 CALL minmax(iip1*jjm, zv,vmin,vllm ) 134 135 ullm = ABS ( ullm ) 136 vllm = ABS ( vllm ) 137 138 zllm = MAX( ullm,vllm ) 139 z1llm = 1./ zllm 140 DO ij = 1, ip1jmp1 141 zu(ij) = zu(ij)* z1llm 142 ENDDO 143 DO ij = 1, ip1jm 144 zv(ij) = zv(ij)* z1llm 145 ENDDO 146 5 CONTINUE 147 148 IF ( ii.EQ.1 ) THEN 149 IF(lstardis) THEN 150 cdivu = 1./zllm 151 ELSE 152 cdivu = zllm **( -1./nitergdiv ) 153 ENDIF 154 ELSE 155 IF(lstardis) THEN 156 crot = 1./ zllm 157 ELSE 158 crot = zllm **( -1./nitergrot ) 159 ENDIF 160 ENDIF 161 162 20 CONTINUE 163 164 c petit test pour les operateurs non star: 165 c ---------------------------------------- 166 167 c IF(.NOT.lstardis) THEN 168 fact = rad*24./REAL(jjm) 169 fact = fact*fact 170 PRINT*,'coef u ', fact/cdivu, 1./cdivu 171 PRINT*,'coef r ', fact/crot , 1./crot 172 PRINT*,'coef h ', fact/cdivh, 1./cdivh 173 c ENDIF 174 175 c----------------------------------------------------------------------- 176 c variation verticale du coefficient de dissipation: 177 c -------------------------------------------------- 178 179 c First step: going from 1 to dissip_fac_mid (in gcm.def) 180 c============ 181 DO l=1,llm 182 zz = 1. - preff/presnivs(l) 183 zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) 184 ENDDO 185 186 write(*,*) 'Dissipation : ' 187 write(*,*) 'Multiplication de la dissipation en altitude :' 188 write(*,*) ' dissip_fac_mid =', dissip_fac_mid 189 190 c Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) 191 c========================== 192 c Utilisation de la fonction tangente hyperbolique pour augmenter 193 c arbitrairement la dissipation et donc la stabilite du modele a 194 c partir d'une certaine altitude. 195 196 c Le facteur multiplicatif de basse atmosphere etant deja pris 197 c en compte, il faut diviser le facteur multiplicatif de haute 198 c atmosphere par celui-ci. 199 200 if (ok_strato) then 201 202 Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) 203 do l=1,llm 204 zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) 205 & *(1-(0.5*(1+tanh(-6./dissip_deltaz* 206 & (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) 207 enddo 208 209 write(*,*) ' dissip_fac_up =', dissip_fac_up 210 write(*,*) 'Transition mid /up: Pupstart,delta =', 211 & dissip_pupstart,'Pa', dissip_deltaz , '(km)' 212 213 endif 214 215 216 PRINT*,'Constantes de temps de la diffusion horizontale' 217 218 tetamin = 1.e+6 219 220 DO l=1,llm 221 tetaudiv(l) = zvert(l)/tetagdiv 222 tetaurot(l) = zvert(l)/tetagrot 223 tetah(l) = zvert(l)/tetatemp 224 225 IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 226 IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 227 IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) 228 ENDDO 229 230 PRINT *,' INIDI tetamin dtvr ',tetamin,dtvr,iperiod 231 idissip = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod 232 PRINT *,' INIDI tetamin idissip ',tetamin,idissip 233 idissip = MAX(iperiod,idissip) 234 dtdiss = idissip * dtvr 235 PRINT *,' INIDI idissip dtdiss ',idissip,dtdiss 236 237 DO l = 1,llm 238 PRINT*,zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), 239 * dtdiss*tetah(l) 240 ENDDO 241 242 c 243 RETURN 244 END 57 ENDDO 58 ENDDO 59 60 idum = -1 61 zh(1) = RAN1(idum)-.5 62 idum = 0 63 DO ij = 2, ip1jmp1 64 zh(ij) = RAN1(idum) -.5 65 ENDDO 66 67 CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1) 68 69 CALL minmax(iip1*jjp1,zh,zhmin,zhmax ) 70 71 IF ( zhmin .GE. zhmax ) THEN 72 write(lunout,*)' Inidissip zh min max ',zhmin,zhmax 73 abort_message='probleme generateur alleatoire dans inidissip' 74 call abort_gcm('inidissip',abort_message,1) 75 ENDIF 76 77 zllm = ABS( zhmax ) 78 DO l = 1,50 79 IF(lstardis) THEN 80 CALL divgrad2(1,zh,deltap,niterh,zh) 81 ELSE 82 CALL divgrad (1,zh,niterh,zh) 83 ENDIF 84 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 92 ENDDO 93 94 IF(lstardis) THEN 95 cdivh = 1./ zllm 96 ELSE 97 cdivh = zllm ** ( -1./niterh ) 98 ENDIF 99 100 ! calcul des valeurs propres de gradiv (ii =1) et nxgrarot(ii=2) 101 ! ----------------------------------------------------------------- 102 write(lunout,*)'inidissip: calcul des valeurs propres' 103 104 DO ii = 1, 2 105 ! 106 DO ij = 1, ip1jmp1 107 zu(ij) = RAN1(idum) -.5 108 ENDDO 109 CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1) 110 DO ij = 1, ip1jm 111 zv(ij) = RAN1(idum) -.5 112 ENDDO 113 CALL filtreg (zv,jjm,1,2,1,.FALSE.,1) 114 115 CALL minmax(iip1*jjp1,zu,umin,ullm ) 116 CALL minmax(iip1*jjm, zv,vmin,vllm ) 117 118 ullm = ABS ( ullm ) 119 vllm = ABS ( vllm ) 120 121 DO l = 1, 50 122 IF(ii.EQ.1) THEN 123 !cccc CALL covcont( 1,zu,zv,zu,zv ) 124 IF(lstardis) THEN 125 CALL gradiv2( 1,zu,zv,nitergdiv,zu,zv ) 126 ELSE 127 CALL gradiv ( 1,zu,zv,nitergdiv,zu,zv ) 128 ENDIF 129 ELSE 130 IF(lstardis) THEN 131 CALL nxgraro2( 1,zu,zv,nitergrot,zu,zv ) 132 ELSE 133 CALL nxgrarot( 1,zu,zv,nitergrot,zu,zv ) 134 ENDIF 135 ENDIF 136 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 151 end DO 152 153 IF ( ii.EQ.1 ) THEN 154 IF(lstardis) THEN 155 cdivu = 1./zllm 156 ELSE 157 cdivu = zllm **( -1./nitergdiv ) 158 ENDIF 159 ELSE 160 IF(lstardis) THEN 161 crot = 1./ zllm 162 ELSE 163 crot = zllm **( -1./nitergrot ) 164 ENDIF 165 ENDIF 166 167 end DO 168 169 ! petit test pour les operateurs non star: 170 ! ---------------------------------------- 171 172 ! IF(.NOT.lstardis) THEN 173 fact = rad*24./REAL(jjm) 174 fact = fact*fact 175 write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu 176 write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot 177 write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh 178 ! ENDIF 179 180 !----------------------------------------------------------------------- 181 ! variation verticale du coefficient de dissipation: 182 ! -------------------------------------------------- 183 184 ! First step: going from 1 to dissip_fac_mid (in gcm.def) 185 !============ 186 DO l=1,llm 187 zz = 1. - preff/presnivs(l) 188 zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz ) 189 ENDDO 190 191 write(lunout,*) 'Dissipation : ' 192 write(lunout,*) 'Multiplication de la dissipation en altitude :' 193 write(lunout,*) ' dissip_fac_mid =', dissip_fac_mid 194 195 ! Second step if ok_strato: from dissip_fac_mid to dissip_fac_up (in gcm.def) 196 !========================== 197 ! Utilisation de la fonction tangente hyperbolique pour augmenter 198 ! arbitrairement la dissipation et donc la stabilite du modele a 199 ! partir d'une certaine altitude. 200 201 ! Le facteur multiplicatif de basse atmosphere etant deja pris 202 ! en compte, il faut diviser le facteur multiplicatif de haute 203 ! atmosphere par celui-ci. 204 205 if (ok_strato) then 206 207 Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta) 208 do l=1,llm 209 zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) & 210 *(1-(0.5*(1+tanh(-6./dissip_deltaz* & 211 (-dissip_hdelta*log(presnivs(l)/Pup)) )))) )) 212 enddo 213 214 write(*,*) ' dissip_fac_up =', dissip_fac_up 215 write(*,*) 'Transition mid /up: Pupstart,delta =', & 216 dissip_pupstart,'Pa', dissip_deltaz , '(km)' 217 218 endif 219 220 221 write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale' 222 223 tetamin = 1.e+6 224 225 DO l=1,llm 226 tetaudiv(l) = zvert(l)/tetagdiv 227 tetaurot(l) = zvert(l)/tetagrot 228 tetah(l) = zvert(l)/tetatemp 229 230 IF( tetamin.GT. (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l) 231 IF( tetamin.GT. (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l) 232 IF( tetamin.GT. (1./ tetah(l)) ) tetamin = 1./ tetah(l) 233 ENDDO 234 235 ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def 236 IF (dissip_period == 0) THEN 237 dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod 238 write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period 239 dissip_period = MAX(iperiod,dissip_period) 240 END IF 241 242 dtdiss = dissip_period * dtvr 243 write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr 244 245 DO l = 1,llm 246 write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), & 247 dtdiss*tetah(l) 248 ENDDO 249 250 END SUBROUTINE inidissip -
trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F
r7 r270 1 1 ! 2 ! $Id: integrd_p.F 1 446 2010-10-22 09:27:25Z emillour$2 ! $Id: integrd_p.F 1550 2011-07-05 09:44:55Z lguez $ 3 3 ! 4 4 SUBROUTINE integrd_p … … 128 128 PRINT*,' Au point ij = ',stop_it, ' , pression sol neg. ' 129 129 & , ps(stop_it) 130 STOP' dans integrd' 130 print *, ' dans integrd' 131 stop 1 131 132 ENDIF 132 133 -
trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F
r130 r270 168 168 169 169 character*80 dynhist_file, dynhistave_file 170 character *20 modname170 character(len=*),parameter :: modname="leapfrog" 171 171 character*80 abort_message 172 172 … … 226 226 endif 227 227 itaufinp1 = itaufin +1 228 modname="leapfrog_p"229 228 230 229 itau = 0 … … 392 391 ! Purely Matsuno time stepping 393 392 IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE. 394 IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE. 393 IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward ) 394 s apdiss = .TRUE. 395 395 IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 396 396 s .and. iflag_phys.EQ.1 ) apphys = .TRUE. … … 398 398 ! Leapfrog/Matsuno time stepping 399 399 IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE. 400 IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE. 400 IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward ) 401 s apdiss = .TRUE. 401 402 IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE. 402 403 END IF … … 650 651 c ------------------------------------------------------------- 651 652 652 IF( forward. OR . leapf ) THEN 653 ! IF( forward. OR . leapf ) THEN 654 IF((.not.forward).OR. leapf ) THEN 655 ! Ehouarn: gather mass fluxes during backward Matsuno or LF step 653 656 cc$OMP PARALLEL DEFAULT(SHARED) 654 657 c -
trunk/LMDZ.COMMON/libf/dyn3dpar/limit_netcdf.F90
r7 r270 1 1 ! 2 ! $Id: limit_netcdf.F90 1 441 2010-10-13 13:06:56Z emillour$2 ! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $ 3 3 !------------------------------------------------------------------------------- 4 4 ! … … 42 42 REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask 43 43 #ifndef CPP_EARTH 44 WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'44 CALL abort_gcm('limit_netcdf','Earth-specific routine, needs Earth physics',1) 45 45 #else 46 46 !------------------------------------------------------------------------------- … … 52 52 #include "indicesol.h" 53 53 54 !--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------55 LOGICAL, PARAMETER :: fracterre=.TRUE.56 57 54 !--- INPUT NETCDF FILES NAMES -------------------------------------------------- 58 55 CHARACTER(LEN=25) :: icefile, sstfile, dumstr 59 56 CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc ', & 60 57 famipsic='amipbc_sic_1x1.nc ', & 61 fclimsst='amipbc_sst_1x1_clim.nc ', &62 fclimsic='amipbc_sic_1x1_clim.nc ', &63 58 fcpldsst='cpl_atm_sst.nc ', & 64 59 fcpldsic='cpl_atm_sic.nc ', & 60 fhistsst='histmth_sst.nc ', & 61 fhistsic='histmth_sic.nc ', & 65 62 frugo ='Rugos.nc ', & 66 63 falbe ='Albedo.nc ' 67 64 CHARACTER(LEN=10) :: varname 68 65 !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------ 69 66 REAL, DIMENSION(klon) :: fi_ice, verif … … 80 77 INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC 81 78 INTEGER :: NF90_FORMAT 82 LOGICAL :: lCPL !--- T: IPCC-IPSL cpl model output files83 79 INTEGER :: ndays !--- Depending on the output calendar 84 80 … … 104 100 105 101 !--- RUGOSITY TREATMENT -------------------------------------------------------- 106 WRITE(lunout,*) 'Traitement de la rugosite' 107 CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 102 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite' 103 varname='RUGOS' 104 CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:)) 108 105 109 106 !--- OCEAN TREATMENT ----------------------------------------------------------- 110 PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.107 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique' 111 108 112 109 ! Input SIC file selection 113 icefile='fake' 114 IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic) 115 IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic) 116 IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic) 117 SELECT CASE(icefile) 118 CASE(famipsic); dumstr='Amip.' 119 CASE(fclimsic); dumstr='Amip climatologique.' 120 CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. 121 CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1) 122 END SELECT 110 ! Open file only to test if available 111 IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 112 icefile=TRIM(famipsic) 113 varname='sicbcs' 114 ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 115 icefile=TRIM(fcpldsic) 116 varname='SIICECOV' 117 ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 118 icefile=TRIM(fhistsic) 119 varname='pourc_sic' 120 ELSE 121 WRITE(lunout,*) 'ERROR! No sea-ice input file was found.' 122 WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ',trim(fhistsic) 123 CALL abort_gcm('limit_netcdf','No sea-ice file was found',1) 124 END IF 123 125 ierr=NF90_CLOSE(nid) 124 WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)125 126 CALL get_2Dfield(icefile, 'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL)126 IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile) 127 128 CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice) 127 129 128 130 ALLOCATE(pctsrf_t(klon,nbsrf,ndays)) 129 131 DO k=1,ndays 130 fi_ice=phy_ice(:,k) 131 WHERE(fi_ice>=1.0 ) fi_ice=1.0 132 WHERE(fi_ice<EPSFRA) fi_ice=0.0 133 IF(fracterre) THEN 134 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil 135 pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice 136 IF(lCPL) THEN ! SIC=pICE*(1-LIC-TER) 137 pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 138 ELSE ! SIC=pICE-LIC 132 fi_ice=phy_ice(:,k) 133 WHERE(fi_ice>=1.0 ) fi_ice=1.0 134 WHERE(fi_ice<EPSFRA) fi_ice=0.0 135 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil 136 pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice 137 IF (icefile==trim(fcpldsic)) THEN ! SIC=pICE*(1-LIC-TER) 138 pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter)) 139 ELSE IF (icefile==trim(fhistsic)) THEN ! SIC=pICE 140 pctsrf_t(:,is_sic,k)=fi_ice(:) 141 ELSE ! icefile==famipsic ! SIC=pICE-LIC 139 142 pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k) 140 141 142 143 END IF 144 WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0. 145 WHERE(1.0-zmasq<EPSFRA) 143 146 pctsrf_t(:,is_sic,k)=0.0 144 147 pctsrf_t(:,is_oce,k)=0.0 145 148 ELSEWHERE 146 149 WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq) 147 pctsrf_t(:,is_sic,k)=1.0-zmasq148 pctsrf_t(:,is_oce,k)=0.0150 pctsrf_t(:,is_sic,k)=1.0-zmasq 151 pctsrf_t(:,is_oce,k)=0.0 149 152 ELSEWHERE 150 pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)151 WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)152 pctsrf_t(:,is_oce,k)=0.0153 pctsrf_t(:,is_sic,k)=1.0-zmasq154 END WHERE153 pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k) 154 WHERE(pctsrf_t(:,is_oce,k)<EPSFRA) 155 pctsrf_t(:,is_oce,k)=0.0 156 pctsrf_t(:,is_sic,k)=1.0-zmasq 157 END WHERE 155 158 END WHERE 156 END WHERE 157 nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) 158 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 159 nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) 160 IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad 161 ELSE 162 pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) 163 WHERE(NINT(pctsrf(:,is_ter))==1) 164 pctsrf_t(:,is_sic,k)=0. 165 pctsrf_t(:,is_oce,k)=0. 166 WHERE(fi_ice>=1.e-5) 167 pctsrf_t(:,is_lic,k)=fi_ice 168 pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k) 169 ELSEWHERE 170 pctsrf_t(:,is_lic,k)=0.0 171 END WHERE 172 ELSEWHERE 173 pctsrf_t(:,is_lic,k) = 0.0 174 WHERE(fi_ice>=1.e-5) 175 pctsrf_t(:,is_ter,k)=0.0 176 pctsrf_t(:,is_sic,k)=fi_ice 177 pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k) 178 ELSEWHERE 179 pctsrf_t(:,is_sic,k)=0.0 180 pctsrf_t(:,is_oce,k)=1.0 181 END WHERE 182 END WHERE 183 verif=sum(pctsrf_t(:,:,k),dim=2) 184 nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5) 185 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 186 END IF 159 END WHERE 160 nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0) 161 IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad 162 nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA) 163 IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad 187 164 END DO 188 165 DEALLOCATE(phy_ice) 189 166 190 167 !--- SST TREATMENT ------------------------------------------------------------- 191 WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.168 IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst' 192 169 193 170 ! Input SST file selection 194 sstfile='fake' 195 IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst) 196 IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst) 197 IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst) 198 SELECT CASE(icefile) 199 CASE(famipsic); dumstr='Amip.' 200 CASE(fclimsic); dumstr='Amip climatologique.' 201 CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE. 202 CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1) 203 END SELECT 171 ! Open file only to test if available 172 IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 173 sstfile=TRIM(famipsst) 174 varname='tosbcs' 175 ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 176 sstfile=TRIM(fcpldsst) 177 varname='SISUTESW' 178 ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN 179 sstfile=TRIM(fhistsst) 180 varname='tsol_oce' 181 ELSE 182 WRITE(lunout,*) 'ERROR! No sst input file was found.' 183 WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst) 184 CALL abort_gcm('limit_netcdf','No sst file was found',1) 185 END IF 204 186 ierr=NF90_CLOSE(nid) 205 WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)206 207 CALL get_2Dfield( trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)187 IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile) 188 189 CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap) 208 190 209 191 !--- ALBEDO TREATMENT ---------------------------------------------------------- 210 WRITE(lunout,*) 'Traitement de l albedo' 211 CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb) 192 IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo' 193 varname='ALBEDO' 194 CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb) 212 195 213 196 !--- REFERENCE GROUND HEAT FLUX TREATMENT -------------------------------------- … … 215 198 216 199 !--- OUTPUT FILE WRITING ------------------------------------------------------- 217 WRITE(lunout,*) 'Ecriture du fichier limit'200 IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut' 218 201 219 202 !--- File creation … … 264 247 ierr=NF90_CLOSE(nid) 265 248 249 IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin' 250 266 251 DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug) 267 252 … … 276 261 !------------------------------------------------------------------------------- 277 262 ! 278 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)263 SUBROUTINE get_2Dfield(fnam, varname, mode, ibar, ndays, champo, flag, mask) 279 264 ! 280 265 !----------------------------------------------------------------------------- … … 304 289 ! Arguments: 305 290 CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name 291 CHARACTER(LEN=10), INTENT(IN) :: varname ! NetCDF variable name 306 292 CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB 307 293 LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels 308 294 INTEGER, INTENT(IN) :: ndays ! current year number of days 309 REAL, POINTER, DIMENSION(:, :) :: champo 295 REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t) 310 296 LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC) 311 297 REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask 312 LOGICAL, OPTIONAL, INTENT(IN) :: lCPL ! Coupled model flag (for ICE)313 298 !------------------------------------------------------------------------------ 314 299 ! Local variables: … … 316 301 INTEGER :: ncid, varid ! NetCDF identifiers 317 302 CHARACTER(LEN=30) :: dnam ! dimension name 318 CHARACTER(LEN=80) :: varname ! NetCDF variable name319 303 !--- dimensions 320 304 INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers … … 331 315 !--- input files 332 316 CHARACTER(LEN=20) :: cal_in ! calendar 317 CHARACTER(LEN=20) :: unit_sic ! attribute unit in sea-ice file 333 318 INTEGER :: ndays_in ! number of days 334 319 !--- misc … … 337 322 CHARACTER(LEN=25) :: title ! for messages 338 323 LOGICAL :: extrp ! flag for extrapolation 324 LOGICAL :: oldice ! flag for old way ice computation 339 325 REAL :: chmin, chmax 340 326 INTEGER ierr 341 327 integer n_extrap ! number of extrapolated points 342 328 logical skip 329 343 330 !------------------------------------------------------------------------------ 344 331 !---Variables depending on keyword 'mode' ------------------------------------- 345 332 NULLIFY(champo) 333 346 334 SELECT CASE(mode) 347 CASE('RUG'); varname='RUGOS';title='Rugosite'348 CASE('SIC'); varname='sicbcs'; title='Sea-ice'349 CASE('SST'); varname='tosbcs'; title='SST'350 CASE('ALB'); varname='ALBEDO'; title='Albedo'335 CASE('RUG'); title='Rugosite' 336 CASE('SIC'); title='Sea-ice' 337 CASE('SST'); title='SST' 338 CASE('ALB'); title='Albedo' 351 339 END SELECT 340 341 352 342 extrp=.FALSE. 343 oldice=.FALSE. 353 344 IF ( PRESENT(flag) ) THEN 354 345 IF ( flag .AND. mode=='SST' ) extrp=.TRUE. 346 IF ( flag .AND. mode=='SIC' ) oldice=.TRUE. 355 347 END IF 356 348 357 349 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ----------------------------- 350 IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam 358 351 ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid); CALL ncerr(ierr, fnam) 359 ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam)352 ierr=NF90_INQ_VARID(ncid, trim(varname), varid); CALL ncerr(ierr, fnam) 360 353 ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam) 354 355 !--- Read unit for sea-ice variable only 356 IF (mode=='SIC') THEN 357 ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic) 358 IF(ierr/=NF90_NOERR) THEN 359 IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value' 360 unit_sic='X' 361 ELSE 362 IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic 363 END IF 364 END IF 361 365 362 366 !--- Longitude … … 365 369 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) 366 370 ierr=NF90_GET_VAR(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam) 367 WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep371 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep 368 372 369 373 !--- Latitude … … 372 376 ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam) 373 377 ierr=NF90_GET_VAR(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam) 374 WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep378 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep 375 379 376 380 !--- Time (variable is not needed - it is rebuilt - but calendar is) … … 385 389 CASE('SIC', 'SST'); cal_in='gregorian' 386 390 END SELECT 387 WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &391 IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' & 388 392 // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.' 389 393 END IF 390 WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &394 IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', & 391 395 cal_in 392 396 397 393 398 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION -------------------- 394 399 !--- Determining input file number of days, depending on calendar … … 398 403 !--- If input records are not monthly, time sampling has to be constant ! 399 404 timeyear=mid_months(anneeref, cal_in, lmdep) 400 IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' & 401 // TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, & 402 ' enregistrements.' 405 IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), & 406 ' ne comportent pas 12, mais ', lmdep, ' enregistrements.' 403 407 404 408 !--- GETTING THE FIELD AND INTERPOLATING IT ---------------------------------- … … 406 410 IF(extrp) ALLOCATE(work(imdep, jmdep)) 407 411 408 WRITE(lunout, *) 409 WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, & 410 ' CHAMPS.' 412 IF (prt_level>5) WRITE(lunout, *) 413 IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.' 411 414 ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam) 412 415 DO l=1, lmdep … … 419 422 work) 420 423 421 IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN 422 IF(l==1) THEN 423 WRITE(lunout, *) & 424 '-------------------------------------------------------------------------' 425 WRITE(lunout, *) & 426 'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$' 427 WRITE(lunout, *) & 428 '-------------------------------------------------------------------------' 424 IF(ibar .AND. .NOT.oldice) THEN 425 IF(l==1 .AND. prt_level>5) THEN 426 WRITE(lunout, *) '-------------------------------------------------------------------------' 427 WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$' 428 WRITE(lunout, *) '-------------------------------------------------------------------------' 429 429 END IF 430 430 IF(mode=='RUG') champ=LOG(champ) … … 453 453 454 454 !--- TIME INTERPOLATION ------------------------------------------------------ 455 WRITE(lunout, *) 456 WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' 457 WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear 458 WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays 455 IF (prt_level>5) THEN 456 WRITE(lunout, *) 457 WRITE(lunout, *)'INTERPOLATION TEMPORELLE.' 458 WRITE(lunout, *)' Vecteur temps en entree: ', timeyear 459 WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays 460 END IF 461 459 462 ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays)) 460 463 skip = .false. … … 471 474 END DO 472 475 if (n_extrap /= 0) then 473 print *,"get_2Dfield pchfe_95: n_extrap = ", n_extrap476 WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap 474 477 end if 475 478 champan(iip1, :, :)=champan(1, :, :) … … 479 482 DO j=1, jjp1 480 483 CALL minmax(iip1, champan(1, j, 10), chmin, chmax) 481 WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j484 IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j 482 485 END DO 483 486 484 487 !--- SPECIAL FILTER FOR SST: SST>271.38 -------------------------------------- 485 488 IF(mode=='SST') THEN 486 WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'489 IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38' 487 490 WHERE(champan<271.38) champan=271.38 488 491 END IF … … 490 493 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ------------------------------------- 491 494 IF(mode=='SIC') THEN 492 WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' 493 IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100. 495 IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0' 496 497 IF (unit_sic=='1') THEN 498 ! Nothing to be done for sea-ice field is already in fraction of 1 499 ! This is the case for sea-ice in file cpl_atm_sic.nc 500 IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1' 501 ELSE 502 ! Convert sea ice from percentage to fraction of 1 503 IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' 504 champan(:, :, :)=champan(:, :, :)/100. 505 END IF 506 494 507 champan(iip1, :, :)=champan(1, :, :) 495 508 WHERE(champan>1.0) champan=1.0 496 509 WHERE(champan<0.0) champan=0.0 497 510 END IF 498 511 499 512 !--- DYNAMICAL TO PHYSICAL GRID ---------------------------------------------- -
trunk/LMDZ.COMMON/makelmdz
r134 r270 389 389 fi 390 390 391 INCLUDE="$INCLUDE -I${NETCDF_INCDIR}"392 LIB="$LIB -L${NETCDF_LIBDIR} -lnetcdf"391 INCLUDE="$INCLUDE ${NETCDF_INCDIR}" 392 LIB="$LIB ${NETCDF_LIBDIR}" 393 393 394 394 ######################################################################## -
trunk/LMDZ.COMMON/makelmdz_fcm
r134 r270 256 256 # LIB="${LIB} -l${LIBPREFIX}cosp" 257 257 fi 258 INCLUDE="$INCLUDE -I${NETCDF_INCDIR}" 259 LIB="$LIB -L${NETCDF_LIBDIR} -lnetcdf" 258 259 INCLUDE="$INCLUDE ${NETCDF_INCDIR}" 260 LIB="$LIB ${NETCDF_LIBDIR}" 260 261 261 262 ########################################################################
Note: See TracChangeset
for help on using the changeset viewer.