Changeset 2193 for trunk/LMDZ.VENUS
- Timestamp:
- Dec 13, 2019, 10:44:42 AM (5 years ago)
- Location:
- trunk/LMDZ.VENUS/libf/phyvenus
- Files:
-
- 3 edited
-
phys_state_var_mod.F90 (modified) (1 diff)
-
physiq_mod.F (modified) (19 diffs)
-
phytrac_chimie.F (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.VENUS/libf/phyvenus/phys_state_var_mod.F90
r2192 r2193 119 119 120 120 !====================================================================== 121 SUBROUTINE phys_state_var_init 121 SUBROUTINE phys_state_var_init(nqmax) 122 122 123 123 IMPLICIT NONE 124 124 #include "dimsoil.h" 125 126 integer :: nqmax 125 127 126 128 ALLOCATE(ftsol(klon)) ! temperature de surface -
trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F
r2135 r2193 90 90 USE mod_phys_lmdz_para, only : is_parallel,jj_nb, 91 91 & is_north_pole_phy, 92 & is_south_pole_phy 92 & is_south_pole_phy, 93 & is_master 93 94 #endif 94 95 IMPLICIT none … … 171 172 172 173 c Variables propres a la physique 173 c 174 INTEGER,save :: itap ! compteur pour la physique174 175 integer,save :: itap ! physics counter 175 176 REAL delp(klon,klev) ! epaisseur d'une couche 176 177 REAL omega(klon,klev) ! vitesse verticale en Pa/s 177 178 178 179 179 INTEGER igwd,idx(klon),itest(klon) 180 c 180 181 181 c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro 182 182 … … 275 275 c 276 276 REAL tmpout(klon,klev) ! K s-1 277 278 INTEGER itaprad279 SAVE itaprad280 277 c 281 278 REAL dist, rmu0(klon), fract(klon) … … 292 289 REAL varerr 293 290 291 ! photochemistry 292 293 integer :: chempas 294 real :: nbapp_chem, zctime 295 296 ! sedimentation 297 298 REAL :: m0_mode1(klev,2),m0_mode2(klev,2) 299 REAL :: m3_mode1(klev,3),m3_mode2 (klev,3) 300 REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2) 301 REAL :: aer_flux(klev) 302 REAL :: d_tr_ssed(klon) 303 c 294 304 c Variables du changement 295 305 c … … 348 358 REAL :: d_tr(klon,klev,nqmax) 349 359 350 c Variables tendance sedimentation351 352 REAL :: m0_mode1(klev,2),m0_mode2(klev,2)353 REAL :: m3_mode1(klev,3),m3_mode2 (klev,3)354 REAL :: d_drop_sed(klev),d_ccn_sed(klev,2),d_liq_sed(klev,2)355 REAL :: aer_flux(klev)356 REAL :: d_tr_sed(klon,klev,nqmax)357 REAL :: d_tr_ssed(klon)358 c359 360 c pour ioipsl 360 361 INTEGER nid_day, nid_mth, nid_ins … … 458 459 . if_ebil) 459 460 460 call phys_state_var_init 461 call phys_state_var_init(nqmax) 461 462 c 462 463 c Initialising Hedin model for upper atm 463 464 c (to be revised when coupled to chemistry) : 464 465 call conc_init 465 c 466 c Initialiser les compteurs: 467 c 468 itap = 0 469 itaprad = 0 466 467 ! initialise physics counter 468 469 itap = 0 470 470 471 471 #ifdef MESOSCALE … … 523 523 ENDIF 524 524 525 radpas = NINT( 525 radpas = NINT(RDAY/pdtphys/nbapp_rad) 526 526 527 527 CALL printflag( ok_journe,ok_instan ) 528 528 529 529 #ifdef CPP_XIOS 530 531 530 write(*,*) "physiq: call initialize_xios_output" 532 531 call initialize_xios_output(rjourvrai,gmtime,pdtphys,RDAY, … … 715 714 endif 716 715 717 c number of microphysical tracers 718 nmicro=0 719 if (ok_cloud .AND. (cl_scheme.eq.1)) nmicro=2 720 if (ok_cloud .AND. (cl_scheme.eq.2)) nmicro=12 716 ! number of microphysical tracers 717 718 nmicro = 0 719 if (ok_cloud .and. (cl_scheme == 1)) nmicro = 2 720 if (ok_cloud .and. (cl_scheme == 2)) nmicro = 12 721 721 722 722 c CAS 1D POUR MICROPHYS Aurelien … … 747 747 748 748 ENDIF ! debut 749 c====================================================================== 750 c====================================================================== 749 751 750 ! ------------------------------------------------------ 752 751 ! Initializations done at every physical timestep: … … 915 914 c 916 915 CALL hgardfou(t_seri,ftsol,'debutphy') 917 c====================================================================918 c919 c Incrementer le compteur de la physique920 c921 itap = itap + 1922 916 923 917 c==================================================================== … … 946 940 ENDIF 947 941 948 c==================================================================== 949 c Calcul des tendances traceurs 950 c==================================================================== 942 ! print fraction of venus day 943 944 if (is_master) then 945 print*, 'gmtime = ', gmtime 946 end if 951 947 952 948 if (iflag_trac.eq.1) then 953 949 !==================================================================== 950 ! Case 1: pseudo-chemistry with relaxation toward fixed profile 951 !==================================================================== 954 952 if (tr_scheme.eq.1) then 955 ! Case 1: pseudo-chemistry with relaxation toward fixed profile956 953 957 954 call phytrac_relax (debut,lafin,nqmax, … … 960 957 961 958 elseif (tr_scheme.eq.2) then 959 !==================================================================== 962 960 ! Case 2: surface emission 963 961 ! For the moment, inspired from Mars version 964 962 ! However, the variable 'source' could be used in physiq 965 963 ! so the call to phytrac_emiss could be to initialise it. 964 !==================================================================== 966 965 967 966 call phytrac_emiss ( (rjourvrai+gmtime)*RDAY, … … 971 970 O tr_seri) 972 971 973 elseif (tr_scheme.eq.3) then ! identical to ok_chem.or.ok_cloud 974 ! Case 3: Full chemistry and/or clouds 975 976 call phytrac_chimie( 977 I debut, 978 I gmtime, 979 I nqmax, 980 I klon, 981 I latitude_deg, 982 I longitude_deg, 983 I nlev, 984 I dtime, 985 I t_seri, 986 I pplay, 987 O tr_seri) 988 989 c CALL WriteField_phy('Pression',pplay,nlev) 990 c CALL WriteField_phy('PressionBnd',paprs,nlev+1) 991 c CALL WriteField_phy('Temp',t_seri,nlev) 992 c IF (ok_cloud) THEN 993 c CALL WriteField_phy('NBRTOT',NBRTOT,nlev) 994 c ENDIF 995 c CALL WriteField_phy('SAl',tr_seri(:,:,i_h2so4liq),nlev) 996 c CALL WriteField_phy('SAg',tr_seri(:,:,i_h2so4),nlev) 997 998 C === SEDIMENTATION === 972 else if (tr_scheme.eq.3) then 973 !==================================================================== 974 ! Case 3: Full chemistry and/or clouds. 975 ! routines are called every "chempas" physical timestep. 976 ! if the physics is called 96000 times per venus day, 977 ! then chempas = 4 978 !==================================================================== 979 980 nbapp_chem = 24000 ! corresponds to 420 seconds 981 chempas = nint(rday/pdtphys/nbapp_chem) 982 zctime = dtime*real(chempas) ! chemical timestep (420s) 983 984 if (mod(itap,chempas) == 0) then ! <------- start of chemistry supercycling 985 986 ! photochemistry and microphysics 987 988 call phytrac_chimie(debut, 989 $ gmtime, 990 $ nqmax, 991 $ klon, 992 $ latitude_deg, 993 $ longitude_deg, 994 $ nlev, 995 $ zctime, 996 $ t_seri, 997 $ pplay, 998 $ tr_seri, 999 $ d_tr_chem, 1000 $ iter) 999 1001 1000 1002 if (ok_sedim) then 1001 1002 c !! Depend on cl_scheme !! 1003 1004 if (cl_scheme.eq.1) then 1005 c ================ 1003 if (cl_scheme == 1) then 1004 1005 ! sedimentation for simplified microphysics 1006 1006 1007 #ifndef MESOSCALE 1007 CALL new_cloud_sedim( 1008 I klon, 1009 I nlev, 1010 I dtime, 1011 I pplay, 1012 I paprs, 1013 I t_seri, 1014 I tr_seri, 1015 O d_tr_sed(:,:,1:2), 1016 O d_tr_ssed, 1017 I nqmax, 1018 O Fsedim) 1019 1020 DO k = 1, klev 1021 DO i = 1, klon 1022 1023 c-------------------- 1024 c Ce test est necessaire pour eviter Xliq=NaN 1025 ! IF (ieee_is_nan(d_tr_sed(i,k,1)).OR. 1026 ! & ieee_is_nan(d_tr_sed(i,k,2))) THEN 1027 IF ((d_tr_sed(i,k,1).ne.d_tr_sed(i,k,1)).OR. 1028 & (d_tr_sed(i,k,2).ne.d_tr_sed(i,k,2))) THEN 1029 PRINT*,'sedim NaN PROBLEM' 1030 PRINT*,'d_tr_sed Nan?',d_tr_sed(i,k,:),'Temp',t_seri(i,k) 1031 PRINT*,'lat-lon',i,'level',k,'dtime',dtime 1032 PRINT*,'F_sed',Fsedim(i,k) 1033 PRINT*,'===============================================' 1034 d_tr_sed(i,k,:)=0. 1035 ENDIF 1036 c-------------------- 1037 1038 tr_seri(i,k,i_h2so4liq) = tr_seri(i,k,i_h2so4liq)+ 1039 & d_tr_sed(i,k,1) 1040 tr_seri(i,k,i_h2oliq) = tr_seri(i,k,i_h2oliq)+ 1041 & d_tr_sed(i,k,2) 1042 d_tr_sed(i,k,:) = d_tr_sed(i,k,:) / dtime 1043 Fsedim(i,k) = Fsedim(i,k) / dtime 1044 1045 ENDDO 1046 ENDDO 1047 1048 Fsedim(:,klev+1) = 0. 1049 1050 elseif (cl_scheme.eq.2) then 1051 c ==================== 1052 1053 d_tr_sed(:,:,:) = 0.D0 1054 1055 DO i=1, klon 1056 1057 c Mode 1 1058 m0_mode1(:,1)=tr_seri(i,:,i_m0_mode1drop) 1059 m0_mode1(:,2)=tr_seri(i,:,i_m0_mode1ccn) 1060 m3_mode1(:,1)=tr_seri(i,:,i_m3_mode1sa) 1061 m3_mode1(:,2)=tr_seri(i,:,i_m3_mode1w) 1062 m3_mode1(:,3)=tr_seri(i,:,i_m3_mode1ccn) 1063 1064 call drop_sedimentation(dtime,klev,m0_mode1,m3_mode1, 1065 & paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:), 1066 & 1, 1067 & d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed) 1008 call new_cloud_sedim(klon, 1009 $ nlev, 1010 $ zctime, 1011 $ pplay, 1012 $ paprs, 1013 $ t_seri, 1014 $ tr_seri, 1015 $ d_tr_sed(:,:,1:2), 1016 $ d_tr_ssed, 1017 $ nqmax, 1018 $ Fsedim) 1019 1020 ! test to avoid nans 1021 1022 do k = 1, klev 1023 do i = 1, klon 1024 if ((d_tr_sed(i,k,1) /= d_tr_sed(i,k,1)) .or. 1025 $ (d_tr_sed(i,k,2) /= d_tr_sed(i,k,2))) then 1026 print*,'sedim NaN PROBLEM' 1027 print*,'d_tr_sed Nan?',d_tr_sed(i,k,:) 1028 print*,'Temp',t_seri(i,k) 1029 print*,'lat-lon',i,'level',k,'zctime',zctime 1030 print*,'F_sed',Fsedim(i,k) 1031 d_tr_sed(i,k,:) = 0. 1032 end if 1033 end do 1034 end do 1035 1036 ! tendency due to sedimentation 1037 1038 d_tr_sed(:,:,:) = d_tr_sed(:,:,:)/zctime 1039 Fsedim(:,1:klev) = Fsedim(:,1:klev)/zctime 1040 Fsedim(:,klev+1) = 0. 1041 1042 else if (cl_scheme == 2) then 1043 1044 ! sedimentation for detailed microphysics 1045 1046 d_tr_sed(:,:,:) = 0. 1047 1048 do i = 1, klon 1049 1050 ! mode 1 1051 1052 m0_mode1(:,1) = tr_seri(i,:,i_m0_mode1drop) 1053 m0_mode1(:,2) = tr_seri(i,:,i_m0_mode1ccn) 1054 m3_mode1(:,1) = tr_seri(i,:,i_m3_mode1sa) 1055 m3_mode1(:,2) = tr_seri(i,:,i_m3_mode1w) 1056 m3_mode1(:,3) = tr_seri(i,:,i_m3_mode1ccn) 1057 1058 call drop_sedimentation(zctime,klev,m0_mode1,m3_mode1, 1059 $ paprs(i,:),zzlev(i,:), 1060 $ zzlay(i,:),t_seri(i,:),1, 1061 $ d_ccn_sed(:,1),d_drop_sed, 1062 $ d_ccn_sed(:,2),d_liq_sed) 1068 1063 1069 1064 d_tr_sed(i,:,i_m0_mode1drop)= d_tr_sed(i,:,i_m0_mode1drop) 1070 &+ d_drop_sed1065 $ + d_drop_sed 1071 1066 d_tr_sed(i,:,i_m0_mode1ccn) = d_tr_sed(i,:,i_m0_mode1ccn) 1072 &+ d_ccn_sed(:,1)1067 $ + d_ccn_sed(:,1) 1073 1068 d_tr_sed(i,:,i_m3_mode1ccn) = d_tr_sed(i,:,i_m3_mode1ccn) 1074 &+ d_ccn_sed(:,2)1069 $ + d_ccn_sed(:,2) 1075 1070 d_tr_sed(i,:,i_m3_mode1sa) = d_tr_sed(i,:,i_m3_mode1sa) 1076 &+ d_liq_sed(:,1)1071 $ + d_liq_sed(:,1) 1077 1072 d_tr_sed(i,:,i_m3_mode1w) = d_tr_sed(i,:,i_m3_mode1w) 1078 & + d_liq_sed(:,2) 1079 1080 c Mode 2 1081 m0_mode2(:,1)=tr_seri(i,:,i_m0_mode2drop) 1082 m0_mode2(:,2)=tr_seri(i,:,i_m0_mode2ccn) 1083 m3_mode2(:,1)=tr_seri(i,:,i_m3_mode2sa) 1084 m3_mode2(:,2)=tr_seri(i,:,i_m3_mode2w) 1085 m3_mode2(:,3)=tr_seri(i,:,i_m3_mode2ccn) 1086 1087 call drop_sedimentation(dtime,klev,m0_mode2,m3_mode2, 1088 & paprs(i,:),zzlev(i,:),zzlay(i,:),t_seri(i,:), 1089 & 2, 1090 & d_ccn_sed(:,1),d_drop_sed,d_ccn_sed(:,2),d_liq_sed) 1073 $ + d_liq_sed(:,2) 1074 1075 ! mode 2 1076 1077 m0_mode2(:,1) = tr_seri(i,:,i_m0_mode2drop) 1078 m0_mode2(:,2) = tr_seri(i,:,i_m0_mode2ccn) 1079 m3_mode2(:,1) = tr_seri(i,:,i_m3_mode2sa) 1080 m3_mode2(:,2) = tr_seri(i,:,i_m3_mode2w) 1081 m3_mode2(:,3) = tr_seri(i,:,i_m3_mode2ccn) 1082 1083 call drop_sedimentation(zctime,klev,m0_mode2,m3_mode2, 1084 $ paprs(i,:),zzlev(i,:), 1085 $ zzlay(i,:),t_seri(i,:),2, 1086 $ d_ccn_sed(:,1),d_drop_sed, 1087 $ d_ccn_sed(:,2),d_liq_sed) 1091 1088 1092 1089 d_tr_sed(i,:,i_m0_mode2drop)= d_tr_sed(i,:,i_m0_mode2drop) 1093 &+ d_drop_sed1090 $ + d_drop_sed 1094 1091 d_tr_sed(i,:,i_m0_mode2ccn) = d_tr_sed(i,:,i_m0_mode2ccn) 1095 &+ d_ccn_sed(:,1)1092 $ + d_ccn_sed(:,1) 1096 1093 d_tr_sed(i,:,i_m3_mode2ccn) = d_tr_sed(i,:,i_m3_mode2ccn) 1097 &+ d_ccn_sed(:,2)1094 $ + d_ccn_sed(:,2) 1098 1095 d_tr_sed(i,:,i_m3_mode2sa) = d_tr_sed(i,:,i_m3_mode2sa) 1099 &+ d_liq_sed(:,1)1096 $ + d_liq_sed(:,1) 1100 1097 d_tr_sed(i,:,i_m3_mode2w) = d_tr_sed(i,:,i_m3_mode2w) 1101 & + d_liq_sed(:,2) 1102 1103 c Aer 1104 call aer_sedimentation(dtime,klev,tr_seri(i,:,i_m0_aer), 1105 & tr_seri(i,:,i_m3_aer),paprs(i,:), 1106 & zzlev(i,:),zzlay(i,:),t_seri(i,:), 1107 & d_tr_sed(i,:,i_m0_aer),d_tr_sed(i,:,i_m3_aer),aer_flux) 1108 1109 END DO 1098 $ + d_liq_sed(:,2) 1099 1100 ! aer 1101 1102 call aer_sedimentation(zctime,klev, 1103 $ tr_seri(i,:,i_m0_aer), 1104 $ tr_seri(i,:,i_m3_aer), 1105 $ paprs(i,:),zzlev(i,:), 1106 $ zzlay(i,:),t_seri(i,:), 1107 $ d_tr_sed(i,:,i_m0_aer), 1108 $ d_tr_sed(i,:,i_m3_aer), 1109 $ aer_flux) 1110 1111 end do 1110 1112 1111 DO iq=nqmax-nmicro+1,nqmax 1112 tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_sed(:,:,iq) 1113 d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq) / dtime 1114 END DO 1113 ! tendency due to sedimentation 1114 1115 do iq = nqmax - nmicro + 1,nqmax 1116 tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_sed(:,:,iq) 1117 d_tr_sed(:,:,iq) = d_tr_sed(:,:,iq)/zctime 1118 end do 1115 1119 #endif 1116 endif 1117 c ==================== 1118 1119 C === FIN SEDIMENTATION === 1120 1121 endif ! ok_sedim 1122 1123 endif ! tr_scheme 1124 endif ! iflag_trac 1125 1126 c 1120 end if ! cl_scheme 1121 end if ! ok_sedim 1122 end if ! mod(itap,chempas) <------- end of chemistry supercycling 1123 1124 ! update tracers (chemistry) 1125 1126 do iq = 1, nqmax - nmicro 1127 tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_chem(:,:,iq)*dtime 1128 end do 1129 1130 ! update tracers (sedimentation) 1131 1132 tr_seri(:,:,i_h2so4liq) = tr_seri(:,:,i_h2so4liq) 1133 $ + d_tr_sed(:,:,1)*dtime 1134 tr_seri(:,:,i_h2oliq) = tr_seri(:,:,i_h2oliq) 1135 $ + d_tr_sed(:,:,2)*dtime 1136 1137 end if ! tr_scheme 1138 end if ! iflag_trac 1139 1127 1140 c==================================================================== 1128 1141 c Appeler la diffusion verticale (programme de couche limite) … … 1337 1350 END IF 1338 1351 1339 1340 c====================================================================1341 c RAYONNEMENT1342 c====================================================================1343 1344 1352 c------------------------------------ 1345 c . Compute radiative tendencies : 1346 c------------------------------------ 1347 c==================================================================== 1348 IF (MOD(itaprad,radpas).EQ.0) THEN 1349 c==================================================================== 1350 1351 dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s) 1352 c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas 1353 1354 1355 c------------------------------------ 1356 c . Compute mean mass, cp and R : 1353 c Compute mean mass, cp and R : 1357 1354 c------------------------------------ 1358 1355 1359 1356 if(callthermos) then 1360 1357 call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax) 1361 1362 1358 endif 1363 1359 1364 1365 cc!!! ADD key callhedin 1366 1367 IF(callnlte.or.callthermos) THEN 1368 call compo_hedin83_mod(pplay,rmu0, 1369 & co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm) 1370 1371 IF(ok_chem) then 1360 c==================================================================== 1361 c RAYONNEMENT 1362 c==================================================================== 1363 if (mod(itap,radpas) == 0) then 1364 1365 dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s) 1366 c==================================================================== 1367 1368 if (callnlte .or. callthermos) then 1369 if (ok_chem) then 1370 1371 ! nlte : use computed chemical species 1372 1372 1373 CC !! GG : Using only mayor species tracers abundances to compute NLTE heating/cooling 1374 1375 CC Conversion [mmr] ---> [vmr] 1376 1377 co2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co2)* 1378 & mmean(1:nlon,1:nlev)/M_tr(i_co2) 1379 covmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_co)* 1380 & mmean(1:nlon,1:nlev)/M_tr(i_co) 1381 ovmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_o)* 1382 & mmean(1:nlon,1:nlev)/M_tr(i_o) 1383 n2vmr_gcm(:,:) = tr_seri(1:nlon,1:nlev,i_n2)* 1384 & mmean(1:nlon,1:nlev)/M_tr(i_n2) 1385 1386 ENDIF 1387 1388 ENDIF 1389 1390 c 1373 co2vmr_gcm(:,:) = tr_seri(:,:,i_co2)*mmean(:,:)/m_tr(i_co2) 1374 covmr_gcm(:,:) = tr_seri(:,:,i_co)*mmean(:,:)/m_tr(i_co) 1375 ovmr_gcm(:,:) = tr_seri(:,:,i_o)*mmean(:,:)/m_tr(i_o) 1376 n2vmr_gcm(:,:) = tr_seri(:,:,i_n2)*mmean(:,:)/m_tr(i_n2) 1377 1378 else 1379 1380 ! nlte : use hedin climatology 1381 1382 call compo_hedin83_mod(pplay,rmu0, 1383 $ co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm) 1384 end if 1385 end if 1386 1391 1387 c NLTE cooling from CO2 emission 1392 1388 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 1406 1402 write(*,*)'I will continue anyway' 1407 1403 endif 1408 1409 1404 endif 1410 1405 … … 1480 1475 c $ covmr_gcm, ovmr_gcm,d_t_euv ) 1481 1476 call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay, 1482 $ rmu0, pdtphys,gmtime,rjourvrai,1477 $ rmu0,dtimerad,gmtime,rjourvrai, 1483 1478 $ tr_seri, d_tr, d_t_euv ) 1484 1479 … … 1493 1488 1494 1489 c==================================================================== 1495 itaprad = 01496 1490 ENDIF ! radpas 1497 1491 c==================================================================== … … 1504 1498 ENDDO 1505 1499 ENDDO 1500 1501 ! increment physics counter 1502 1503 itap = itap + 1 1506 1504 1507 1505 ! CONDUCTION and MOLECULAR VISCOSITY … … 1978 1976 CALL send_xios_field("fluxec",flux_ec) 1979 1977 1980 c When using tracers 1981 IF (iflag_trac.eq.1) THEN 1982 c photochemical compounds !!!outputs in [vmr] 1983 DO iq=1,nqmax-nmicro 1984 CALL send_xios_field(tname(iq),qx(:,:,iq)*mmean(:,:)/M_tr(iq)) 1985 ENDDO 1986 1987 IF ((tr_scheme.eq.3).and.(cl_scheme.eq.1)) THEN 1988 c liquids !!!outputs in [vmr] 1989 CALL send_xios_field(tname(i_h2oliq), 1990 . qx(:,:,i_h2oliq)*mmean(:,:)/M_tr(i_h2oliq)) 1991 CALL send_xios_field(tname(i_h2so4liq), 1992 . qx(:,:,i_h2so4liq)*mmean(:,:)/M_tr(i_h2so4liq)) 1993 if (ok_sedim) CALL send_xios_field("Fsedim",Fsedim(:,1:klev)) 1994 ENDIF 1995 ENDIF 1978 ! when using tracers 1979 1980 if (iflag_trac == 1) then 1981 1982 ! tracers in gas phase, volume mixing ratio 1983 1984 do iq = 1,nqmax - nmicro 1985 call send_xios_field(tname(iq), 1986 $ qx(:,:,iq)*mmean(:,:)/m_tr(iq)) 1987 end do 1988 1989 ! tracers in liquid phase, volume mixing ratio 1990 1991 if ((tr_scheme == 3) .and. (cl_scheme == 1)) THEN 1992 call send_xios_field(tname(i_h2oliq), 1993 $ qx(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq)) 1994 call send_xios_field(tname(i_h2so4liq), 1995 $ qx(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq)) 1996 if (ok_sedim) then 1997 call send_xios_field("Fsedim",fsedim(:,1:klev)) 1998 end if 1999 end if 2000 2001 ! chemical iterations 2002 2003 call send_xios_field("iter",real(iter)) 2004 2005 end if 1996 2006 1997 2007 IF (callthermos .and. ok_chem) THEN -
trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F
r2188 r2193 3 3 $ gmtime, 4 4 $ nqmax, 5 $ n _lon,5 $ nlon, 6 6 $ lat, 7 7 $ lon, 8 $ n _lev,8 $ nlev, 9 9 $ pdtphys, 10 10 $ temp, 11 11 $ pplay, 12 $ trac )13 14 ! $ iter) temporary 12 $ trac, 13 $ d_tr_chem, 14 $ iter) 15 15 16 16 use chemparam_mod … … 26 26 !=================================================================== 27 27 28 integer :: n_lon, n_lev ! number of gridpoints and levels 29 integer :: nqmax ! number of tracers 30 31 real :: gmtime 32 real :: pdtphys ! physics timestep (s) 33 real, dimension(n_lon,n_lev) :: temp ! temperature (k) 34 real, dimension(n_lon,n_lev) :: pplay ! pressure (pa) 35 36 logical :: debutphy ! first call flag 28 integer :: nlon, nlev ! number of gridpoints and levels 29 integer :: nqmax ! number of tracers 30 31 real :: gmtime ! day fraction 32 real :: pdtphys ! phytrac_chimie timestep (s) 33 real, dimension(nlon,nlev) :: temp ! temperature (k) 34 real, dimension(nlon,nlev) :: pplay ! pressure (pa) 35 real, dimension(nlon,nlev,nqmax) :: trac ! tracer mass mixing ratio 36 37 logical :: debutphy ! first call flag 37 38 38 39 !=================================================================== … … 40 41 !=================================================================== 41 42 42 integer, dimension(n_lon,n_lev) :: iter ! chemical iterations 43 44 !=================================================================== 45 ! input/output 46 !=================================================================== 47 48 real, dimension(n_lon,n_lev,nqmax) :: trac ! tracer mass mixing ratio 43 real, dimension(nlon,nlev,nqmax) :: d_tr_chem ! chemical tendency for each tracer 44 integer, dimension(nlon,nlev) :: iter ! chemical iterations 49 45 50 46 !=================================================================== … … 58 54 integer :: ilon, ilev 59 55 60 real lat(n_lon), lat_local(n_lon) 61 real lon(n_lon), lon_local(n_lon) 62 63 real, dimension(n_lon,n_lev) :: mrtwv, mrtsa ! total water and total sulfuric acid 64 real, dimension(n_lon,n_lev) :: mrwv, mrsa ! gas-phase water and gas-phase sulfuric acid 65 real, dimension(n_lon,n_lev) :: trac_sum 56 real lat(nlon), lat_local(nlon) 57 real lon(nlon), lon_local(nlon) 58 59 real, dimension(nlon,nlev) :: mrtwv, mrtsa ! total water and total sulfuric acid 60 real, dimension(nlon,nlev) :: mrwv, mrsa ! gas-phase water and gas-phase sulfuric acid 61 real, dimension(nlon,nlev) :: trac_sum 62 real, dimension(nlon,nlev,nqmax) :: ztrac ! local tracer mixing ratio 66 63 67 64 !=================================================================== … … 120 117 ! case of detailed microphysics without chemistry 121 118 !------------------------------------------------------------------- 122 if (.not. ok_chem .and. ok_cloud .and. cl_scheme ==2) then119 if (.not. ok_chem .and. ok_cloud .and. cl_scheme == 2) then 123 120 124 121 ! convert mass to volume mixing ratio 125 122 126 123 do iq = 1,nqmax - nmicro 127 trac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)124 ztrac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq) 128 125 end do 129 126 130 127 ! initialise microphysics 131 128 132 call vapors4muphy_ini(n _lon,n_lev,trac)129 call vapors4muphy_ini(nlon,nlev,ztrac) 133 130 134 131 ! convert volume to mass mixing ratio 135 132 136 133 do iq = 1,nqmax - nmicro 137 trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)134 trac(:,:,iq) = ztrac(:,:,iq)*m_tr(iq)/mmean(:,:) 138 135 end do 139 136 … … 147 144 148 145 do iq = 1,nqmax - nmicro 149 trac(:,:,iq) = max(trac(:,:,iq)*mmean(:,:)/m_tr(iq),1.e-30)146 ztrac(:,:,iq) = max(trac(:,:,iq)*mmean(:,:)/m_tr(iq), 1.e-30) 150 147 end do 151 148 … … 158 155 ! convert mass to volume mixing ratio : liquid phase 159 156 160 trac(:,:,i_h2so4liq) = max(trac(:,:,i_h2so4liq)161 $ *mmean(:,:)/m_tr(i_h2so4liq),1.e-30)162 trac(:,:,i_h2oliq) = max(trac(:,:,i_h2oliq)163 $ *mmean(:,:)/m_tr(i_h2oliq),1.e-30)157 ztrac(:,:,i_h2so4liq) = max(trac(:,:,i_h2so4liq) 158 $ *mmean(:,:)/m_tr(i_h2so4liq), 1.e-30) 159 ztrac(:,:,i_h2oliq) = max(trac(:,:,i_h2oliq) 160 $ *mmean(:,:)/m_tr(i_h2oliq), 1.e-30) 164 161 165 162 ! total water and sulfuric acid (gas + liquid) 166 163 167 mrtwv(:,:) = trac(:,:,i_h2o) +trac(:,:,i_h2oliq)168 mrtsa(:,:) = trac(:,:,i_h2so4) +trac(:,:,i_h2so4liq)164 mrtwv(:,:) = ztrac(:,:,i_h2o) + ztrac(:,:,i_h2oliq) 165 mrtsa(:,:) = ztrac(:,:,i_h2so4) + ztrac(:,:,i_h2so4liq) 169 166 170 167 ! all water and sulfuric acid is put in the gas-phase … … 175 172 ! call microphysics 176 173 177 call new_cloud_venus(n _lev, n_lon, temp, pplay,174 call new_cloud_venus(nlev, nlon, temp, pplay, 178 175 $ mrtwv, mrtsa, mrwv, mrsa) 179 176 180 177 ! update water vapour and sulfuric acid 181 178 182 trac(:,:,i_h2o) = mrwv(:,:)183 trac(:,:,i_h2oliq) = mrtwv(:,:) -trac(:,:,i_h2o)179 ztrac(:,:,i_h2o) = mrwv(:,:) 180 ztrac(:,:,i_h2oliq) = mrtwv(:,:) - ztrac(:,:,i_h2o) 184 181 185 trac(:,:,i_h2so4) = mrsa(:,:)186 trac(:,:,i_h2so4liq) = mrtsa(:,:) -trac(:,:,i_h2so4)182 ztrac(:,:,i_h2so4) = mrsa(:,:) 183 ztrac(:,:,i_h2so4liq) = mrtsa(:,:) - ztrac(:,:,i_h2so4) 187 184 188 185 end if ! simplified scheme … … 194 191 if (ok_cloud .and. cl_scheme == 2) then 195 192 196 c Boucle sur grille (n_lon) et niveaux (n_lev) 197 DO ilon=1, n_lon 198 DO ilev=1, n_lev 199 200 if (temp(ilon,ilev).lt.500.) then 201 CALL MAD_MUPHY(pdtphys, ! Timestep 202 & temp(ilon,ilev),pplay(ilon,ilev), ! Temperature and pressure 203 & trac(ilon,ilev,i_h2o),trac(ilon,ilev,i_h2so4), ! Mixing ratio of SA and W 204 & trac(ilon,ilev,i_m0_aer),trac(ilon,ilev,i_m3_aer), ! Moments of aerosols 205 & trac(ilon,ilev,i_m0_mode1drop),trac(ilon,ilev,i_m0_mode1ccn), ! Moments of mode 1 206 & trac(ilon,ilev,i_m3_mode1sa),trac(ilon,ilev,i_m3_mode1w), ! Moments of mode 1 207 & trac(ilon,ilev,i_m3_mode1ccn), ! Moments of mode 1 208 & trac(ilon,ilev,i_m0_mode2drop),trac(ilon,ilev,i_m0_mode2ccn), ! Moments of mode 2 209 & trac(ilon,ilev,i_m3_mode2sa),trac(ilon,ilev,i_m3_mode2w), ! Moments of mode 2 210 & trac(ilon,ilev,i_m3_mode2ccn)) ! Moments of mode 2 211 else 212 trac(ilon,ilev,i_m0_aer)=0. 213 trac(ilon,ilev,i_m3_aer)=0. 214 trac(ilon,ilev,i_m0_mode1drop)=0. 215 trac(ilon,ilev,i_m0_mode1ccn)=0. 216 trac(ilon,ilev,i_m3_mode1sa)=0. 217 trac(ilon,ilev,i_m3_mode1w)=0. 218 trac(ilon,ilev,i_m3_mode1ccn)=0. 219 trac(ilon,ilev,i_m0_mode2drop)=0. 220 trac(ilon,ilev,i_m0_mode2ccn)=0. 221 trac(ilon,ilev,i_m3_mode2sa)=0. 222 trac(ilon,ilev,i_m3_mode2w)=0. 223 trac(ilon,ilev,i_m3_mode2ccn)=0. 224 endif 225 ENDDO 226 ENDDO 193 do ilon = 1,nlon 194 do ilev = 1, nlev 195 if (temp(ilon,ilev) < 500.) then 196 call mad_muphy(pdtphys, ! timestep 197 $ temp(ilon,ilev),pplay(ilon,ilev), ! temperature and pressure 198 $ ztrac(ilon,ilev,i_h2o), 199 $ ztrac(ilon,ilev,i_h2so4), 200 $ ztrac(ilon,ilev,i_m0_aer), 201 $ ztrac(ilon,ilev,i_m3_aer), 202 $ ztrac(ilon,ilev,i_m0_mode1drop), 203 $ ztrac(ilon,ilev,i_m0_mode1ccn), 204 $ ztrac(ilon,ilev,i_m3_mode1sa), 205 $ ztrac(ilon,ilev,i_m3_mode1w), 206 $ ztrac(ilon,ilev,i_m3_mode1ccn), 207 $ ztrac(ilon,ilev,i_m0_mode2drop), 208 $ ztrac(ilon,ilev,i_m0_mode2ccn), 209 $ ztrac(ilon,ilev,i_m3_mode2sa), 210 $ ztrac(ilon,ilev,i_m3_mode2w), 211 $ ztrac(ilon,ilev,i_m3_mode2ccn)) 212 else 213 ztrac(ilon,ilev,i_m0_aer) = 0. 214 ztrac(ilon,ilev,i_m3_aer) = 0. 215 ztrac(ilon,ilev,i_m0_mode1drop) = 0. 216 ztrac(ilon,ilev,i_m0_mode1ccn) = 0. 217 ztrac(ilon,ilev,i_m3_mode1sa) = 0. 218 ztrac(ilon,ilev,i_m3_mode1w) = 0. 219 ztrac(ilon,ilev,i_m3_mode1ccn) = 0. 220 ztrac(ilon,ilev,i_m0_mode2drop) = 0. 221 ztrac(ilon,ilev,i_m0_mode2ccn) = 0. 222 ztrac(ilon,ilev,i_m3_mode2sa) = 0. 223 ztrac(ilon,ilev,i_m3_mode2w) = 0. 224 ztrac(ilon,ilev,i_m3_mode2ccn) = 0. 225 end if 226 end do 227 end do 227 228 228 229 end if ! detailed scheme … … 238 239 lat_local(:) = lat(:)*rpi/180. 239 240 240 do ilon = 1,n _lon241 do ilon = 1,nlon 241 242 242 243 ! solar zenith angle … … 246 247 $ *sin(lon_local(ilon))*sin(lon_sun))*180./rpi 247 248 248 call photochemistry_venus(n _lev, n_lon, pdtphys,249 call photochemistry_venus(nlev, nlon, pdtphys, 249 250 $ pplay(ilon,:)/100., 250 251 $ temp(ilon,:), 251 $ trac(ilon,:,:),252 $ ztrac(ilon,:,:), 252 253 $ mmean(ilon,:), 253 254 $ sza_local, nqmax, iter(ilon,:)) … … 258 259 259 260 !=================================================================== 260 ! co nvert volume to mass mixing ratio261 ! compute tendencies 261 262 !=================================================================== 262 263 … … 264 265 265 266 do iq = 1,nqmax - nmicro 266 trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:) 267 ztrac(:,:,iq) = ztrac(:,:,iq)*m_tr(iq)/mmean(:,:) 268 d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys 267 269 end do 268 270 … … 270 272 271 273 if (ok_cloud .and. cl_scheme == 1) then 272 trac(:,:,i_h2so4liq) = trac(:,:,i_h2so4liq)*m_tr(i_h2so4liq) 273 & /mmean(:,:) 274 trac(:,:,i_h2oliq) = trac(:,:,i_h2oliq)*m_tr(i_h2oliq) 275 & /mmean(:,:) 274 ztrac(:,:,i_h2so4liq) = ztrac(:,:,i_h2so4liq)*m_tr(i_h2so4liq) 275 $ /mmean(:,:) 276 ztrac(:,:,i_h2oliq) = ztrac(:,:,i_h2oliq)*m_tr(i_h2oliq) 277 $ /mmean(:,:) 278 d_tr_chem(:,:,i_h2so4liq) = (ztrac(:,:,i_h2so4liq) 279 $ - trac(:,:,i_h2so4liq))/pdtphys 280 d_tr_chem(:,:,i_h2oliq) = (ztrac(:,:,i_h2oliq) 281 $ - trac(:,:,i_h2oliq))/pdtphys 276 282 end if 277 283 278 284 end subroutine phytrac_chimie
Note: See TracChangeset
for help on using the changeset viewer.
