Changeset 3693 for LMDZ6/trunk/libf
- Timestamp:
- Jun 3, 2020, 6:15:16 PM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90
r3686 r3693 190 190 real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm) 191 191 real :: d_u_nudge(llm),d_v_nudge(llm) 192 real :: du_adv(llm),dv_adv(llm)193 real :: d u_age(llm),dv_age(llm)192 ! real :: d_u_adv(llm),d_v_adv(llm) 193 real :: d_u_age(llm),d_v_age(llm) 194 194 real :: alpha 195 195 real :: ttt … … 273 273 d_u_nudge(:)=0. 274 274 d_v_nudge(:)=0. 275 d u_adv(:)=0.276 d v_adv(:)=0.277 d u_age(:)=0.278 d v_age(:)=0.275 d_u_adv(:)=0. 276 d_v_adv(:)=0. 277 d_u_age(:)=0. 278 d_v_age(:)=0. 279 279 280 280 ! Initialization of Common turb_forcing … … 840 840 841 841 #include "1D_interp_cases.h" 842 ! Vertical advection843 ! call lstendH(llm,nqtot,omega,d_t_vert_adv,d_q_vert_adv,q,temp,u,v,play)844 ! print*,'B d_t_adv ',d_t_adv(1:20)*86400845 ! print*,'B d_t_vert_adv ',d_t_vert_adv(1:20)*86400846 ! print*,'B dt omega ',omega847 842 848 843 !--------------------------------------------------------------------- … … 872 867 teta=temp*(pzero/play)**rkappa 873 868 do l=2,llm-1 869 ! vertical tendencies computed as d X / d t = -W d X / d z 874 870 d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1)) 875 871 d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1)) 876 d_t_vert_adv(l)=-(w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)))/(pzero/play(l))**rkappa 872 ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa 873 d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa 877 874 d_q_vert_adv(l,1)=-w_adv(l)*(q(l+1,1)-q(l-1,1))/(z_adv(l+1)-z_adv(l-1)) 878 875 enddo 876 d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:) 877 d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:) 879 878 d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:) 880 879 d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1) … … 938 937 fcoriolis=2.*sin(rpi*xlat/180.)*romega 939 938 940 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!941 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!942 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!943 ! if (forcing_radconv .or. forcing_fire) then944 ! fcoriolis=0.0945 ! dt_cooling=0.0946 ! d_t_adv=0.0947 ! d_q_adv=0.0948 ! endif949 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!950 951 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!952 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!953 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!954 ! if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice &955 ! & .or.forcing_amma .or. forcing_type.eq.101) then956 ! fcoriolis=0.0 ; ug=0. ; vg=0.957 ! endif958 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!959 960 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!961 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!962 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!963 ! if(forcing_rico) then964 ! dt_cooling=0.965 ! endif966 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!967 968 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!969 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!970 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!971 !CRio:Attention modif sp??cifique cas de Caroline972 ! if (forcing_type==-1) then973 ! fcoriolis=0.974 !975 !on calcule dt_cooling976 ! do l=1,llm977 ! if (play(l).ge.20000.) then978 ! dt_cooling(l)=-1.5/86400.979 ! elseif ((play(l).ge.10000.).and.((play(l).lt.20000.))) then980 ! dt_cooling(l)=-1.5/86400.*(play(l)-10000.)/(10000.)-1./86400.*(20000.-play(l))/10000.*(temp(l)-200.)981 ! else982 ! dt_cooling(l)=-1.*(temp(l)-200.)/86400.983 ! endif984 ! enddo985 !986 ! endif987 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!988 989 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!990 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!991 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!992 ! if (forcing_sandu) then993 ! ug(1:llm)=u_mod(1:llm)994 ! vg(1:llm)=v_mod(1:llm)995 ! endif996 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!997 998 939 IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', & 999 940 fcoriolis, xlat,mxcalc 1000 941 1001 ! print *,'u-ug=',u-ug1002 1003 ! !!!!!!!!!!!!!!!!!!!!!!!1004 ! Geostrophic wind 1005 ! Le calcul ci dessous est insuffisamment precis 1006 ! du_age(1:mxcalc)=fcoriolis*(v(1:mxcalc)-vg(1:mxcalc)) 1007 ! dv_age(1:mxcalc)=-fcoriolis*(u(1:mxcalc)-ug(1:mxcalc)) 1008 !!!!!!!!!!!!!!!!!!!!!!!! 942 !--------------------------------------------------------------------- 943 ! Geostrophic forcing 944 !--------------------------------------------------------------------- 945 946 IF ( forc_geo == 0 ) THEN 947 d_u_age(1:mxcalc)=0. 948 d_v_age(1:mxcalc)=0. 949 ELSE 1009 950 sfdt = sin(0.5*fcoriolis*timestep) 1010 951 cfdt = cos(0.5*fcoriolis*timestep) 1011 952 ! print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep 1012 953 ! 1013 d u_age(1:mxcalc)= -2.*sfdt/timestep* &954 d_u_age(1:mxcalc)= -2.*sfdt/timestep* & 1014 955 & (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) - & 1015 956 & cfdt*(v(1:mxcalc)-vg(1:mxcalc)) ) 1016 957 !! : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc)) 1017 958 ! 1018 d v_age(1:mxcalc)= -2.*sfdt/timestep* &959 d_v_age(1:mxcalc)= -2.*sfdt/timestep* & 1019 960 & (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) + & 1020 961 & sfdt*(v(1:mxcalc)-vg(1:mxcalc)) ) 1021 962 !! : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc)) 1022 ! 1023 !!!!!!!!!!!!!!!!!!!!!!!! 963 ENDIF 964 ! 965 !--------------------------------------------------------------------- 1024 966 ! Nudging 1025 ! !!!!!!!!!!!!!!!!!!!!!!!967 !--------------------------------------------------------------------- 1026 968 d_t_nudge(:) = 0. 1027 969 d_u_nudge(:) = 0. … … 1039 981 ENDDO 1040 982 983 !--------------------------------------------------------------------- 984 ! Optional outputs 985 !--------------------------------------------------------------------- 1041 986 #ifdef OUTPUT_PHYS_SCM 1042 987 CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv) … … 1056 1001 #endif 1057 1002 1058 !1059 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1060 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!1061 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1062 ! if (forcing_fire) THEN1063 ! print*,'Enlever cette section rapidement'1064 ! stop1065 !1066 !1067 !!let ww=if ( alt le 1100 ) then alt*-0.00001 else 01068 !!let wt=if ( alt le 1100 ) then min( -3.75e-5 , -7.5e-8*alt) else 01069 !!let wq=if ( alt le 1100 ) then max( 1.5e-8 , 3e-11*alt) else 01070 ! d_t_adv=0.1071 ! d_q_adv=0.1072 ! teta=temp*(pzero/play)**rkappa1073 ! d_t_adv=0.1074 ! d_q_adv=0.1075 ! do l=2,llm-11076 ! if (zlay(l)<=1100) then1077 ! wwww=-0.00001*zlay(l)1078 ! d_t_adv(l)=-wwww*(teta(l)-teta(l+1))/(zlay(l)-zlay(l+1)) /(pzero/play(l))**rkappa1079 ! d_q_adv(l,1:2)=-wwww*(q(l,1:2)-q(l+1,1:2))/(zlay(l)-zlay(l+1))1080 ! d_t_adv(l)=d_t_adv(l)+min(-3.75e-5 , -7.5e-8*zlay(l))1081 ! d_q_adv(l,1)=d_q_adv(l,1)+max( 1.5e-8 , 3e-11*zlay(l))1082 ! endif1083 ! enddo1084 !1085 ! endif1086 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1087 1088 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1089 ! call writefield_phy('dv_age' ,dv_age,llm)1090 ! call writefield_phy('du_age' ,du_age,llm)1091 ! call writefield_phy('du_phys' ,du_phys,llm)1092 ! call writefield_phy('u_tend' ,u,llm)1093 ! call writefield_phy('u_g' ,ug,llm)1094 !1095 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!1096 !! Increment state variables1097 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!1098 1003 IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added 1099 1004 1100 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1101 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!1102 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1103 ! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h1104 ! au dessus de 700hpa, on relaxe vers les profils initiaux1105 ! if (forcing_sandu .OR. forcing_astex) then1106 !#include "1D_nudge_sandu_astex.h"1107 ! else1108 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1109 1005 u(1:mxcalc)=u(1:mxcalc) + timestep*( & 1110 1006 & du_phys(1:mxcalc) & 1111 & +d u_age(1:mxcalc)+du_adv(1:mxcalc) &1007 & +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc) & 1112 1008 & +d_u_nudge(1:mxcalc) ) 1113 1009 v(1:mxcalc)=v(1:mxcalc) + timestep*( & 1114 1010 & dv_phys(1:mxcalc) & 1115 & +d v_age(1:mxcalc)+dv_adv(1:mxcalc) &1011 & +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc) & 1116 1012 & +d_v_nudge(1:mxcalc) ) 1117 1013 q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*( & … … 1125 1021 & temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) 1126 1022 print* ,'dv_phys=',dv_phys 1127 print* ,'d v_age=',dv_age1128 print* ,'d v_adv=',dv_adv1023 print* ,'d_v_age=',d_v_age 1024 print* ,'d_v_adv=',d_v_adv 1129 1025 print* ,'d_v_nudge=',d_v_nudge 1130 1026 print*, v
Note: See TracChangeset
for help on using the changeset viewer.