Ignore:
Timestamp:
Jun 3, 2020, 6:15:16 PM (4 years ago)
Author:
fhourdin
Message:

Corrections on the standard SCM format.
Fredho

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/dyn1d/scm.F90

    r3686 r3693  
    190190      real :: dt_cooling(llm),d_t_adv(llm),d_t_nudge(llm)
    191191      real :: d_u_nudge(llm),d_v_nudge(llm)
    192       real :: du_adv(llm),dv_adv(llm)
    193       real :: du_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)
    194194      real :: alpha
    195195      real :: ttt
    … …  
    273273      d_u_nudge(:)=0.
    274274      d_v_nudge(:)=0.
    275       du_adv(:)=0.
    276       dv_adv(:)=0.
    277       du_age(:)=0.
    278       dv_age(:)=0.
     275      d_u_adv(:)=0.
     276      d_v_adv(:)=0.
     277      d_u_age(:)=0.
     278      d_v_age(:)=0.
    279279     
    280280! Initialization of Common turb_forcing
    … …  
    840840
    841841#include "1D_interp_cases.h"
    842    ! Vertical advection
    843 !  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)*86400
    845 !  print*,'B d_t_vert_adv ',d_t_vert_adv(1:20)*86400
    846 !  print*,'B dt omega ',omega
    847842
    848843!---------------------------------------------------------------------
    … …  
    872867      teta=temp*(pzero/play)**rkappa
    873868      do l=2,llm-1
     869        ! vertical tendencies computed as d X / d t = -W  d X / d z
    874870        d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l-1))/(z_adv(l+1)-z_adv(l-1))
    875871        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
    877874        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))
    878875      enddo
     876      d_u_adv(:)=d_u_adv(:)+d_u_vert_adv(:)
     877      d_v_adv(:)=d_v_adv(:)+d_v_vert_adv(:)
    879878      d_t_adv(:)=d_t_adv(:)+d_t_vert_adv(:)
    880879      d_q_adv(:,1)=d_q_adv(:,1)+d_q_vert_adv(:,1)
    … …  
    938937       fcoriolis=2.*sin(rpi*xlat/180.)*romega
    939938
    940 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    941 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    942 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    943 !       if (forcing_radconv .or. forcing_fire) then
    944 !         fcoriolis=0.0
    945 !         dt_cooling=0.0
    946 !         d_t_adv=0.0
    947 !         d_q_adv=0.0
    948 !       endif
    949 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    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) then
    956 !        fcoriolis=0.0 ; ug=0. ; vg=0.
    957 !      endif
    958 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    959 
    960 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    961 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    962 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    963 !      if(forcing_rico) then
    964 !         dt_cooling=0.
    965 !      endif
    966 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    967 
    968 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    969 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    970 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    971 !CRio:Attention modif sp??cifique cas de Caroline
    972 !     if (forcing_type==-1) then
    973 !        fcoriolis=0.
    974 !       
    975 !on calcule dt_cooling
    976 !       do l=1,llm
    977 !       if (play(l).ge.20000.) then
    978 !           dt_cooling(l)=-1.5/86400.
    979 !       elseif ((play(l).ge.10000.).and.((play(l).lt.20000.))) then
    980 !           dt_cooling(l)=-1.5/86400.*(play(l)-10000.)/(10000.)-1./86400.*(20000.-play(l))/10000.*(temp(l)-200.)
    981 !       else
    982 !           dt_cooling(l)=-1.*(temp(l)-200.)/86400.
    983 !       endif
    984 !       enddo
    985 !
    986 !     endif     
    987 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    988 
    989 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    990 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    991 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    992 !     if (forcing_sandu) then
    993 !        ug(1:llm)=u_mod(1:llm)
    994 !        vg(1:llm)=v_mod(1:llm)
    995 !     endif
    996 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    997 
    998939      IF (prt_level >= 5) print*, 'fcoriolis, xlat,mxcalc ', &
    999940                                   fcoriolis, xlat,mxcalc
    1000941
    1001 !       print *,'u-ug=',u-ug
    1002 
    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
    1009950       sfdt = sin(0.5*fcoriolis*timestep)
    1010951       cfdt = cos(0.5*fcoriolis*timestep)
    1011952!       print *,'fcoriolis,sfdt,cfdt,timestep',fcoriolis,sfdt,cfdt,timestep
    1012953!
    1013         du_age(1:mxcalc)= -2.*sfdt/timestep*                                &
     954        d_u_age(1:mxcalc)= -2.*sfdt/timestep*                                &
    1014955     &          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -                          &
    1015956     &           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
    1016957!!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
    1017958!
    1018        dv_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
     959       d_v_age(1:mxcalc)= -2.*sfdt/timestep*                                 &
    1019960     &          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +                           &
    1020961     &           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
    1021962!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
    1022 !
    1023 !!!!!!!!!!!!!!!!!!!!!!!!
     963      ENDIF
     964!
     965!---------------------------------------------------------------------
    1024966!  Nudging
    1025 !!!!!!!!!!!!!!!!!!!!!!!!
     967!---------------------------------------------------------------------
    1026968      d_t_nudge(:) = 0.
    1027969      d_u_nudge(:) = 0.
    … …  
    1039981      ENDDO
    1040982
     983!---------------------------------------------------------------------
     984!  Optional outputs
     985!---------------------------------------------------------------------
    1041986#ifdef OUTPUT_PHYS_SCM
    1042987      CALL iophys_ecrit('w_adv',klev,'w_adv','K/day',w_adv)
    … …  
    10561001#endif
    10571002
    1058 !
    1059 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1060 !! CONSERVE EN ATTENDANT QUE LE CAS EN QUESTION FONCTIONNE EN STD !!
    1061 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1062 !       if (forcing_fire) THEN
    1063 !               print*,'Enlever cette section rapidement'
    1064 !               stop
    1065 !               
    1066 !
    1067 !!let ww=if ( alt le 1100 ) then alt*-0.00001 else 0
    1068 !!let wt=if ( alt le 1100 ) then min( -3.75e-5 , -7.5e-8*alt)  else 0
    1069 !!let wq=if ( alt le 1100 ) then max( 1.5e-8 , 3e-11*alt)  else 0
    1070 !           d_t_adv=0.
    1071 !           d_q_adv=0.
    1072 !           teta=temp*(pzero/play)**rkappa
    1073 !           d_t_adv=0.
    1074 !           d_q_adv=0.
    1075 !           do l=2,llm-1
    1076 !              if (zlay(l)<=1100) then
    1077 !                  wwww=-0.00001*zlay(l)
    1078 !                  d_t_adv(l)=-wwww*(teta(l)-teta(l+1))/(zlay(l)-zlay(l+1)) /(pzero/play(l))**rkappa
    1079 !                  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 !              endif
    1083 !           enddo
    1084 !
    1085 !        endif
    1086 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    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 variables
    1097 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    10981003    IF (flag_inhib_forcing == 0) then ! if tendency of forcings should be added
    10991004
    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.h
    1104 ! au dessus de 700hpa, on relaxe vers les profils initiaux
    1105 !     if (forcing_sandu .OR. forcing_astex) then
    1106 !#include "1D_nudge_sandu_astex.h"
    1107 !      else
    1108 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    11091005        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
    11101006     &              du_phys(1:mxcalc)                                       &
    1111      &             +du_age(1:mxcalc)+du_adv(1:mxcalc)                       &
     1007     &             +d_u_age(1:mxcalc)+d_u_adv(1:mxcalc)                       &
    11121008     &             +d_u_nudge(1:mxcalc) )           
    11131009        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
    11141010     &              dv_phys(1:mxcalc)                                       &
    1115      &             +dv_age(1:mxcalc)+dv_adv(1:mxcalc)                       &
     1011     &             +d_v_age(1:mxcalc)+d_v_adv(1:mxcalc)                       &
    11161012     &             +d_v_nudge(1:mxcalc) )
    11171013        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &
    … …  
    11251021     &              temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1)
    11261022           print* ,'dv_phys=',dv_phys
    1127            print* ,'dv_age=',dv_age
    1128            print* ,'dv_adv=',dv_adv
     1023           print* ,'d_v_age=',d_v_age
     1024           print* ,'d_v_adv=',d_v_adv
    11291025           print* ,'d_v_nudge=',d_v_nudge
    11301026           print*, v
Note: See TracChangeset for help on using the changeset viewer.