- Timestamp:
- Jul 24, 2024, 2:54:37 PM (2 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_old_lmdz1d.F90
r5112 r5116 8 8 SUBROUTINE old_lmdz1d 9 9 10 USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar, getin10 USE ioipsl, ONLY: ju2ymds, ymds2ju, ioconf_calendar, getin 11 11 USE phys_state_var_mod, ONLY: phys_state_var_init, phys_state_var_end, & 12 12 clwcon, detr_therm, & … … 68 68 69 69 integer, parameter :: ngrid = 1 70 real:: zcufi = 1.71 real:: zcvfi = 1.72 73 !- real:: nat_surf70 REAL :: zcufi = 1. 71 REAL :: zcvfi = 1. 72 73 !- REAL :: nat_surf 74 74 !- logical :: ok_flux_surf 75 !- real:: fsens76 !- real:: flat77 !- real:: tsurf78 !- real:: rugos79 !- real:: qsol(1:2)80 !- real:: qsurf81 !- real:: psurf82 !- real:: zsurf83 !- real:: albedo75 !- REAL :: fsens 76 !- REAL :: flat 77 !- REAL :: tsurf 78 !- REAL :: rugos 79 !- REAL :: qsol(1:2) 80 !- REAL :: qsurf 81 !- REAL :: psurf 82 !- REAL :: zsurf 83 !- REAL :: albedo 84 84 !- 85 !- real:: time = 0.86 !- real:: time_ini87 !- real:: xlat88 !- real:: xlon89 !- real:: wtsurf90 !- real:: wqsurf91 !- real:: restart_runoff92 !- real:: xagesno93 !- real:: qsolinp94 !- real:: zpicinp85 !- REAL :: time = 0. 86 !- REAL :: time_ini 87 !- REAL :: xlat 88 !- REAL :: xlon 89 !- REAL :: wtsurf 90 !- REAL :: wqsurf 91 !- REAL :: restart_runoff 92 !- REAL :: xagesno 93 !- REAL :: qsolinp 94 !- REAL :: zpicinp 95 95 !- 96 real:: fnday97 real:: day, daytime98 real:: day199 real:: heure100 integer:: jour101 integer:: mois102 integer:: an96 REAL :: fnday 97 REAL :: day, daytime 98 REAL :: day1 99 REAL :: heure 100 INTEGER :: jour 101 INTEGER :: mois 102 INTEGER :: an 103 103 104 104 !--------------------------------------------------------------------- … … 106 106 !--------------------------------------------------------------------- 107 107 108 integer:: kmax = llm108 INTEGER :: kmax = llm 109 109 integer llm700, nq1, nq2 110 110 INTEGER, PARAMETER :: nlev_max = 1000, nqmx = 1000 … … 117 117 real qprof(nlev_max, nqmx) 118 118 119 ! integer:: forcing_type119 ! INTEGER :: forcing_type 120 120 logical :: forcing_les = .FALSE. 121 121 logical :: forcing_armcu = .FALSE. … … 136 136 logical :: forcing_case2 = .FALSE. 137 137 logical :: forcing_SCM = .FALSE. 138 integer:: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file138 INTEGER :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file 139 139 ! (cf read_tsurf1d.F) 140 140 … … 159 159 ! Declarations related to nudging 160 160 !--------------------------------------------------------------------- 161 integer:: nudge_max161 INTEGER :: nudge_max 162 162 parameter (nudge_max = 9) 163 integer:: inudge_RHT = 1164 integer:: inudge_UV = 2163 INTEGER :: inudge_RHT = 1 164 INTEGER :: inudge_UV = 2 165 165 logical :: nudge(nudge_max) 166 real:: t_targ(llm)167 real:: rh_targ(llm)168 real:: u_targ(llm)169 real:: v_targ(llm)166 REAL :: t_targ(llm) 167 REAL :: rh_targ(llm) 168 REAL :: u_targ(llm) 169 REAL :: v_targ(llm) 170 170 171 171 !--------------------------------------------------------------------- 172 172 ! Declarations related to vertical discretization: 173 173 !--------------------------------------------------------------------- 174 real:: pzero = 1.e5175 real:: play (llm), zlay (llm), sig_s(llm), plev(llm + 1)176 real:: playd(llm), zlayd(llm), ap_amma(llm + 1), bp_amma(llm + 1)174 REAL :: pzero = 1.e5 175 REAL :: play (llm), zlay (llm), sig_s(llm), plev(llm + 1) 176 REAL :: playd(llm), zlayd(llm), ap_amma(llm + 1), bp_amma(llm + 1) 177 177 178 178 !--------------------------------------------------------------------- … … 180 180 !--------------------------------------------------------------------- 181 181 182 real:: phi(llm)183 real:: teta(llm), tetal(llm), temp(llm), u(llm), v(llm), w(llm)182 REAL :: phi(llm) 183 REAL :: teta(llm), tetal(llm), temp(llm), u(llm), v(llm), w(llm) 184 184 REAL rot(1, llm) ! relative vorticity, in s-1 185 real:: rlat_rad(1), rlon_rad(1)186 real:: omega(llm + 1), omega2(llm), rho(llm + 1)187 real:: ug(llm), vg(llm), fcoriolis188 real:: sfdt, cfdt189 real:: du_phys(llm), dv_phys(llm), dt_phys(llm)190 real:: dt_dyn(llm)191 real:: dt_cooling(llm), d_t_adv(llm), d_t_nudge(llm)192 real:: d_u_nudge(llm), d_v_nudge(llm)193 real:: du_adv(llm), dv_adv(llm)194 real:: du_age(llm), dv_age(llm)195 real:: alpha196 real:: ttt185 REAL :: rlat_rad(1), rlon_rad(1) 186 REAL :: omega(llm + 1), omega2(llm), rho(llm + 1) 187 REAL :: ug(llm), vg(llm), fcoriolis 188 REAL :: sfdt, cfdt 189 REAL :: du_phys(llm), dv_phys(llm), dt_phys(llm) 190 REAL :: dt_dyn(llm) 191 REAL :: dt_cooling(llm), d_t_adv(llm), d_t_nudge(llm) 192 REAL :: d_u_nudge(llm), d_v_nudge(llm) 193 REAL :: du_adv(llm), dv_adv(llm) 194 REAL :: du_age(llm), dv_age(llm) 195 REAL :: alpha 196 REAL :: ttt 197 197 198 198 REAL, ALLOCATABLE, DIMENSION(:, :) :: q … … 206 206 ! Initialization of surface variables 207 207 !--------------------------------------------------------------------- 208 real:: run_off_lic_0(1)209 real:: fder(1), snsrf(1, nbsrf), qsurfsrf(1, nbsrf)210 real:: tsoil(1, nsoilmx, nbsrf)211 ! real:: agesno(1,nbsrf)208 REAL :: run_off_lic_0(1) 209 REAL :: fder(1), snsrf(1, nbsrf), qsurfsrf(1, nbsrf) 210 REAL :: tsoil(1, nsoilmx, nbsrf) 211 ! REAL :: agesno(1,nbsrf) 212 212 213 213 !--------------------------------------------------------------------- … … 215 215 !--------------------------------------------------------------------- 216 216 logical :: ok_writedem = .TRUE. 217 real:: sollw_in = 0.218 real:: solsw_in = 0.217 REAL :: sollw_in = 0. 218 REAL :: solsw_in = 0. 219 219 220 220 !--------------------------------------------------------------------- … … 223 223 logical :: firstcall = .TRUE. 224 224 logical :: lastcall = .FALSE. 225 real:: phis(1) = 0.0226 real:: dpsrf(1)225 REAL :: phis(1) = 0.0 226 REAL :: dpsrf(1) 227 227 228 228 !--------------------------------------------------------------------- … … 243 243 ! Fichiers et d'autres variables 244 244 !--------------------------------------------------------------------- 245 integer:: k, l, i, it = 1, mxcalc246 integer:: nsrf245 INTEGER :: k, l, i, it = 1, mxcalc 246 INTEGER :: nsrf 247 247 integer jcode 248 248 INTEGER read_climoz 249 249 250 integer:: it_end ! iteration number of the last call250 INTEGER :: it_end ! iteration number of the last call 251 251 !Al1 252 252 integer ecrit_slab_oc !1=ecrit,-1=lit,0=no file … … 293 293 close(1) 294 294 !Al1 295 write(*, *) 'lmdz1d.def lu => unicol.def'295 WRITE(*, *) 'lmdz1d.def lu => unicol.def' 296 296 297 297 ! forcing_type defines the way the SCM is forced: … … 484 484 485 485 ! calend = 'earth_365d' 486 if (calend == 'earth_360d') then486 if (calend == 'earth_360d') THEN 487 487 CALL ioconf_calendar('360_day') 488 write(*, *)'CALENDRIER CHOISI: Terrestre a 360 jours/an'489 else if (calend == 'earth_365d') then488 WRITE(*, *)'CALENDRIER CHOISI: Terrestre a 360 jours/an' 489 else if (calend == 'earth_365d') THEN 490 490 CALL ioconf_calendar('noleap') 491 write(*, *)'CALENDRIER CHOISI: Terrestre a 365 jours/an'492 else if (calend == 'earth_366d') then491 WRITE(*, *)'CALENDRIER CHOISI: Terrestre a 365 jours/an' 492 else if (calend == 'earth_366d') THEN 493 493 CALL ioconf_calendar('all_leap') 494 write(*, *)'CALENDRIER CHOISI: Terrestre bissextile'495 else if (calend == 'gregorian') then494 WRITE(*, *)'CALENDRIER CHOISI: Terrestre bissextile' 495 else if (calend == 'gregorian') THEN 496 496 stop 'gregorian calend should not be used by normal user' 497 497 CALL ioconf_calendar('gregorian') ! not to be used by normal users 498 write(*, *)'CALENDRIER CHOISI: Gregorien'498 WRITE(*, *)'CALENDRIER CHOISI: Gregorien' 499 499 else 500 500 write (*, *) 'ERROR : unknown calendar ', calend … … 509 509 ! Le numero du jour est dans "day". L heure est traitee separement. 510 510 ! La date complete est dans "daytime" (l'unite est le jour). 511 if (nday>0) then511 if (nday>0) THEN 512 512 fnday = nday 513 513 else … … 627 627 CALL phys_state_var_init(read_climoz) 628 628 629 if (ngrid/=klon) then629 if (ngrid/=klon) THEN 630 630 PRINT*, 'stop in inifis' 631 631 PRINT*, 'Probleme de dimensions :' … … 640 640 !! surf_Planck = 0. 641 641 !! surf_Conv = 0. 642 !! write(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv642 !! WRITE(*,*) 'Gateaux-dif Planck,Conv:',surf_Planck,surf_Conv 643 643 !!! a mettre dans le lmdz1d.def ou autre 644 644 !! … … 692 692 IF (forcing_type == 59) THEN 693 693 ! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m 694 write(*, *) '***********************'694 WRITE(*, *) '***********************' 695 695 do l = 1, llm 696 write(*, *) 'l,play(l),presnivs(l): ', l, play(l), presnivs(l)697 if (trouve_700 .and. play(l)<=70000) then696 WRITE(*, *) 'l,play(l),presnivs(l): ', l, play(l), presnivs(l) 697 if (trouve_700 .and. play(l)<=70000) THEN 698 698 llm700 = l 699 699 print *, 'llm700,play=', llm700, play(l) / 100. … … 701 701 endif 702 702 enddo 703 write(*, *) '***********************'703 WRITE(*, *) '***********************' 704 704 ENDIF 705 705 … … 710 710 INCLUDE "old_1D_read_forc_cases.h" 711 711 712 IF (forcing_GCM2SCM) then712 IF (forcing_GCM2SCM) THEN 713 713 write (*, *) 'forcing_GCM2SCM not yet implemented' 714 714 stop 'in initialization' … … 773 773 ! On le met juste avant pour avoir acces a tous les champs 774 774 775 IF (ok_writedem) then 776 775 IF (ok_writedem) THEN 777 776 !-------------------------------------------------------------------------- 778 777 ! pbl_surface_init (called here) and pbl_surface_final (called by phyredem) … … 813 812 PRINT*, 'avant phyredem' 814 813 pctsrf(1, :) = 0. 815 if (nat_surf==0.) then814 if (nat_surf==0.) THEN 816 815 pctsrf(1, is_oce) = 1. 817 816 pctsrf(1, is_ter) = 0. 818 817 pctsrf(1, is_lic) = 0. 819 818 pctsrf(1, is_sic) = 0. 820 else if (nat_surf == 1) then819 else if (nat_surf == 1) THEN 821 820 pctsrf(1, is_oce) = 0. 822 821 pctsrf(1, is_ter) = 1. 823 822 pctsrf(1, is_lic) = 0. 824 823 pctsrf(1, is_sic) = 0. 825 else if (nat_surf == 2) then824 else if (nat_surf == 2) THEN 826 825 pctsrf(1, is_oce) = 0. 827 826 pctsrf(1, is_ter) = 0. 828 827 pctsrf(1, is_lic) = 1. 829 828 pctsrf(1, is_sic) = 0. 830 else if (nat_surf == 3) then829 else if (nat_surf == 3) THEN 831 830 pctsrf(1, is_oce) = 0. 832 831 pctsrf(1, is_ter) = 0. … … 857 856 pbl_tke(:, 2, :) = 1.e-2 858 857 PRINT *, ' pbl_tke dans lmdz1d ' 859 if (prt_level >= 5) then858 if (prt_level >= 5) THEN 860 859 DO nsrf = 1, 4 861 860 PRINT *, 'pbl_tke(1,:,', nsrf, ') ', pbl_tke(1, :, nsrf) … … 937 936 CALL getin('iflag_physiq', iflag_physiq) 938 937 939 if (.not.restart) then938 if (.not.restart) THEN 940 939 iflag_pbl = 5 941 940 CALL phyredem ("startphy.nc") … … 978 977 CALL phys_state_var_end 979 978 !Al1 980 IF (restart) then979 IF (restart) THEN 981 980 PRINT*, 'CALL to restart dyn 1d' 982 981 Call dyn1deta0("start1dyn.nc", plev, play, phi, phis, presnivs, & … … 1007 1006 END IF 1008 1007 !Al1 ================ end restart ================================= 1009 IF (ecrit_slab_oc==1) then1008 IF (ecrit_slab_oc==1) THEN 1010 1009 open(97, file = 'div_slab.dat', STATUS = 'UNKNOWN') 1011 elseif (ecrit_slab_oc==0) then1010 elseif (ecrit_slab_oc==0) THEN 1012 1011 open(97, file = 'div_slab.dat', STATUS = 'OLD') 1013 1012 END IF … … 1016 1015 ! Initialize target profile for RHT nudging if needed 1017 1016 !--------------------------------------------------------------------- 1018 IF (nudge(inudge_RHT)) then1017 IF (nudge(inudge_RHT)) THEN 1019 1018 CALL nudge_RHT_init(plev, play, temp, q(:, 1), t_targ, rh_targ) 1020 1019 END IF 1021 IF (nudge(inudge_UV)) then1020 IF (nudge(inudge_UV)) THEN 1022 1021 CALL nudge_UV_init(plev, play, u, v, u_targ, v_targ) 1023 1022 END IF … … 1034 1033 DO while(it<=it_end) 1035 1034 1036 if (prt_level>=1) then1035 if (prt_level>=1) THEN 1037 1036 PRINT*, 'XXXXXXXXXXXXXXXXXXX ITAP,day,time=', & 1038 1037 it, day, time, it_end, day_step … … 1058 1057 INCLUDE "old_1D_interp_cases.h" 1059 1058 1060 IF (forcing_GCM2SCM) then1059 IF (forcing_GCM2SCM) THEN 1061 1060 write (*, *) 'forcing_GCM2SCM not yet implemented' 1062 1061 stop 'in time loop' … … 1076 1075 ! Listing output for debug prt_level>=1 1077 1076 !--------------------------------------------------------------------- 1078 IF (prt_level>=1) then1077 IF (prt_level>=1) THEN 1079 1078 print *, ' avant physiq : -------- day time ', day, time 1080 write(*, *) 'firstcall,lastcall,phis', &1079 WRITE(*, *) 'firstcall,lastcall,phis', & 1081 1080 firstcall, lastcall, phis 1082 1081 end if 1083 IF (prt_level>=5) then1084 write(*, '(a10,2a4,4a13)') 'BEFOR1 IT=', 'it', 'l', &1082 IF (prt_level>=5) THEN 1083 WRITE(*, '(a10,2a4,4a13)') 'BEFOR1 IT=', 'it', 'l', & 1085 1084 'presniv', 'plev', 'play', 'phi' 1086 write(*, '(a10,2i4,4f13.2)') ('BEFOR1 IT= ', it, l, &1085 WRITE(*, '(a10,2i4,4f13.2)') ('BEFOR1 IT= ', it, l, & 1087 1086 presnivs(l), plev(l), play(l), phi(l), l = 1, llm) 1088 write(*, '(a11,2a4,a11,6a8)') 'BEFOR2', 'it', 'l', &1087 WRITE(*, '(a11,2a4,a11,6a8)') 'BEFOR2', 'it', 'l', & 1089 1088 'presniv', 'u', 'v', 'temp', 'q1', 'q2', 'omega2' 1090 write(*, '(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ', it, l, &1089 WRITE(*, '(a11,2i4,f11.2,5f8.2,e10.2)') ('BEFOR2 IT= ', it, l, & 1091 1090 presnivs(l), u(l), v(l), temp(l), q(l, 1), q(l, 2), omega2(l), l = 1, llm) 1092 1091 END IF … … 1105 1104 ! Listing output for debug 1106 1105 !--------------------------------------------------------------------- 1107 IF (prt_level>=5) then1108 write(*, '(a11,2a4,4a13)') 'AFTER1 IT=', 'it', 'l', &1106 IF (prt_level>=5) THEN 1107 WRITE(*, '(a11,2a4,4a13)') 'AFTER1 IT=', 'it', 'l', & 1109 1108 'presniv', 'plev', 'play', 'phi' 1110 write(*, '(a11,2i4,4f13.2)') ('AFTER1 it= ', it, l, &1109 WRITE(*, '(a11,2i4,4f13.2)') ('AFTER1 it= ', it, l, & 1111 1110 presnivs(l), plev(l), play(l), phi(l), l = 1, llm) 1112 write(*, '(a11,2a4,a11,6a8)') 'AFTER2', 'it', 'l', &1111 WRITE(*, '(a11,2a4,a11,6a8)') 'AFTER2', 'it', 'l', & 1113 1112 'presniv', 'u', 'v', 'temp', 'q1', 'q2', 'omega2' 1114 write(*, '(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ', it, l, &1113 WRITE(*, '(a11,2i4,f11.2,5f8.2,e10.2)') ('AFTER2 it= ', it, l, & 1115 1114 presnivs(l), u(l), v(l), temp(l), q(l, 1), q(l, 2), omega2(l), l = 1, llm) 1116 write(*, '(a11,2a4,a11,5a8)') 'AFTER3', 'it', 'l', &1115 WRITE(*, '(a11,2a4,a11,5a8)') 'AFTER3', 'it', 'l', & 1117 1116 'presniv', 'du_phys', 'dv_phys', 'dt_phys', 'dq1', 'dq2' 1118 write(*, '(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ', it, l, &1117 WRITE(*, '(a11,2i4,f11.2,5f8.2)') ('AFTER3 it= ', it, l, & 1119 1118 presnivs(l), 86400 * du_phys(l), 86400 * dv_phys(l), & 1120 1119 86400 * dt_phys(l), 86400 * dq(l, 1), dq(l, 2), l = 1, llm) 1121 write(*, *) 'dpsrf', dpsrf1120 WRITE(*, *) 'dpsrf', dpsrf 1122 1121 END IF 1123 1122 !--------------------------------------------------------------------- … … 1126 1125 1127 1126 fcoriolis = 2. * sin(rpi * xlat / 180.) * romega 1128 IF (forcing_radconv .or. forcing_fire) then1127 IF (forcing_radconv .or. forcing_fire) THEN 1129 1128 fcoriolis = 0.0 1130 1129 dt_cooling = 0.0 … … 1135 1134 1136 1135 IF (forcing_toga .or. forcing_GCSSold .or. forcing_twpice & 1137 .or.forcing_amma .or. forcing_type==101) then1136 .or.forcing_amma .or. forcing_type==101) THEN 1138 1137 fcoriolis = 0.0 ; ug = 0. ; vg = 0. 1139 1138 END IF 1140 1139 1141 IF(forcing_rico) then1140 IF(forcing_rico) THEN 1142 1141 dt_cooling = 0. 1143 1142 END IF 1144 1143 1145 1144 !CRio:Attention modif sp??cifique cas de Caroline 1146 IF (forcing_type==-1) then1145 IF (forcing_type==-1) THEN 1147 1146 fcoriolis = 0. 1148 1147 !Nudging … … 1150 1149 !on calcule dt_cooling 1151 1150 do l = 1, llm 1152 if (play(l)>=20000.) then1151 if (play(l)>=20000.) THEN 1153 1152 dt_cooling(l) = -1.5 / 86400. 1154 elseif ((play(l)>=10000.).and.((play(l)<20000.))) then1153 elseif ((play(l)>=10000.).and.((play(l)<20000.))) THEN 1155 1154 dt_cooling(l) = -1.5 / 86400. * (play(l) - 10000.) / (10000.) - 1. / 86400. * (20000. - play(l)) / 10000. * (temp(l) - 200.) 1156 1155 else … … 1161 1160 END IF 1162 1161 !RC 1163 IF (forcing_sandu) then1162 IF (forcing_sandu) THEN 1164 1163 ug(1:llm) = u_mod(1:llm) 1165 1164 vg(1:llm) = v_mod(1:llm) … … 1198 1197 d_u_nudge(:) = 0. 1199 1198 d_v_nudge(:) = 0. 1200 IF (nudge(inudge_RHT)) then1199 IF (nudge(inudge_RHT)) THEN 1201 1200 CALL nudge_RHT(timestep, plev, play, t_targ, rh_targ, temp, q(:, 1), & 1202 1201 d_t_nudge, d_q_nudge(:, 1)) 1203 1202 END IF 1204 IF (nudge(inudge_UV)) then1203 IF (nudge(inudge_UV)) THEN 1205 1204 CALL nudge_UV(timestep, plev, play, u_targ, v_targ, u, v, & 1206 1205 d_u_nudge, d_v_nudge) … … 1218 1217 d_q_adv = 0. 1219 1218 do l = 2, llm - 1 1220 if (zlay(l)<=1100) then1219 if (zlay(l)<=1100) THEN 1221 1220 wwww = -0.00001 * zlay(l) 1222 1221 d_t_adv(l) = -wwww * (teta(l) - teta(l + 1)) / (zlay(l) - zlay(l + 1)) / (pzero / play(l))**rkappa … … 1243 1242 ! pour les cas sandu et astex, on reclacule u,v,q,temp et teta dans 1D_nudge_sandu_astex.h 1244 1243 ! au dessus de 700hpa, on relaxe vers les profils initiaux 1245 if (forcing_sandu .OR. forcing_astex) then1244 if (forcing_sandu .OR. forcing_astex) THEN 1246 1245 INCLUDE "1D_nudge_sandu_astex.h" 1247 1246 else … … 1259 1258 + d_q_nudge(1:mxcalc, :)) 1260 1259 1261 if (prt_level>=3) then1260 if (prt_level>=3) THEN 1262 1261 print *, & 1263 1262 'physiq-> temp(1),dt_phys(1),d_t_adv(1),dt_cooling(1) ', & … … 1316 1315 ! Air temperature : 1317 1316 !--------------------------------------------------------------------- 1318 IF (lastcall) then1317 IF (lastcall) THEN 1319 1318 PRINT*, 'Pas de temps final ', it 1320 1319 CALL ju2ymds(daytime, an, mois, jour, heure)
Note: See TracChangeset
for help on using the changeset viewer.
