Changeset 3022 for trunk/LMDZ.COMMON/libf/evolution
- Timestamp:
- Aug 1, 2023, 4:11:55 PM (17 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 3 edited
-
orbit_param_criterion_mod.F90 (modified) (6 diffs)
-
pem.F90 (modified) (1 diff)
-
recomp_orb_param_mod.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/orbit_param_criterion_mod.F90
r2986 r3022 3 3 !======================================================================= 4 4 ! 5 ! Purpose: Compute the maximum n imber of iteration the PEM can do based6 ! on thestopping criterion given by the orbital parameters5 ! Purpose: Compute the maximum number of iteration for PEM based on the 6 ! stopping criterion given by the orbital parameters 7 7 ! 8 ! Author: RV 8 ! Author: RV, JBC 9 9 !======================================================================= 10 10 … … 48 48 49 49 real :: Year ! Year of the simulation 50 real :: timeperi_ls! Year of the simulation51 integer nlask, ilask!,last_ilask ! Loop variable50 real :: Lsp ! Year of the simulation 51 integer nlask, ilask, iylask ! Loop variables 52 52 parameter (nlask = 20001) ! Number of line in Laskar file 53 53 … … 56 56 real max_obl,max_ex,max_lsp ! Maximum value of orbit param given the acceptable percentage 57 57 real min_obl,min_ex,min_lsp ! Maximum value of orbit param given the acceptable percentage 58 59 real max_obl_iter,max_ex_iter,max_lsp_iter ! Maximum year iteration before reaching an unacceptable value 60 61 logical obl_not_found, ex_not_found,lsp_not_found ! Loop variable (first call) 62 logical max_lsp_modulo,min_lsp_modulo ! Variable to check if we are close to the 360 degree modulo for lsp parameter 63 58 real max_obl_iter,max_ex_iter,max_lsp_iter ! Maximum year iteration before reaching an unacceptable value 64 59 real xa,xb,ya,yb,yc 65 60 … … 68 63 ! ********************************************************************** 69 64 70 Year =year_bp_ini+year_PEM! We compute the current year71 timeperi_ls=timeperi*360/2/pi! We convert in degree72 73 call ini_lask_param_mod(nlask) 74 75 print 76 print 77 print 78 print 79 print 80 print *, "orbit_param_criterion, Current lsp=", timeperi_ls65 Year = year_bp_ini + year_PEM ! We compute the current year 66 Lsp = 360. - timeperi*360./(2.*pi) ! We convert in degree 67 68 call ini_lask_param_mod(nlask) ! Allocation of variables 69 70 print*, "orbit_param_criterion, Year in pem.def=", year_bp_ini 71 print*, "orbit_param_criterion, Year in the startpem.nc =", year_PEM 72 print*, "orbit_param_criterion, Current year=", Year 73 print*, "orbit_param_criterion, Current obl=", obliquit 74 print*, "orbit_param_criterion, Current ex=", e_elips 75 print*, "orbit_param_criterion, Current lsp=", Lsp 81 76 82 77 ! We read the file … … 94 89 close(73) 95 90 96 print *, "Coresponding line in the ob_ex_lsp.asc file=", last_ilask 97 98 ! 5% max change case 99 100 max_change_obl=0.05 101 max_change_ex=0.05 102 max_change_lsp=0.05 103 104 max_obl=obliquit*(1+max_change_obl) 105 min_obl=obliquit*(1-max_change_obl) 106 107 max_ex=e_elips*(1+max_change_ex) 108 min_ex=e_elips*(1-max_change_ex) 109 110 max_lsp=timeperi_ls*(1+max_change_lsp) 111 min_lsp=timeperi_ls*(1-max_change_lsp) 112 113 !End of 5% max change case 91 print*, "Coresponding line in the ob_ex_lsp.asc file=", last_ilask 114 92 115 93 !Constant max change case 116 94 117 max_change_obl =0.5118 max_change_ex =0.1119 max_change_lsp =40.95 max_change_obl = 0.5 96 max_change_ex = 0.1 97 max_change_lsp = 40. 120 98 121 99 CALL getin('max_change_obl', max_change_obl) 100 max_obl = obliquit + max_change_obl 101 min_obl = obliquit - max_change_obl 102 print*, "Maximum obliquity accepted=", max_obl 103 print*, "Minimum obliquity accepted=", min_obl 122 104 123 105 CALL getin('max_change_ex', max_change_ex) 106 max_ex = e_elips + max_change_ex 107 min_ex = e_elips - max_change_ex 108 print*, "Maximum excentricity accepted=", max_ex 109 print*, "Minimum excentricity accepted=", min_ex 124 110 125 111 CALL getin('max_change_lsp', max_change_lsp) 126 127 max_obl=obliquit+max_change_obl 128 min_obl=obliquit-max_change_obl 129 130 max_ex=e_elips+max_change_ex 131 min_ex=e_elips-max_change_ex 132 133 max_lsp=timeperi_ls+max_change_lsp 134 min_lsp=timeperi_ls-max_change_lsp 135 136 max_lsp_modulo=.false. 137 min_lsp_modulo=.false. 138 if(max_lsp.gt.360) max_lsp_modulo=.true. 139 if(min_lsp.lt.0) min_lsp_modulo=.true. 112 max_lsp = Lsp + max_change_lsp 113 min_lsp = Lsp - max_change_lsp 114 print*, "Maximum lsp accepted=", max_lsp 115 print*, "Minimum lsp accepted=", min_lsp 140 116 141 117 !End Constant max change case … … 143 119 ! If we do not want some orb parameter to change, they should not be a stopping criterion, 144 120 ! So the number of iteration corresponding is set to maximum 145 if(.not.var_obl) then 146 obl_not_found=.FALSE. 147 else 148 obl_not_found=.TRUE. 121 max_obl_iter = 999999 122 max_ex_iter = 999999 123 max_lsp_iter = 999999 124 125 ! Tendency of the orbital parameter for the considered year gives the limitation between min and max 126 do ilask = last_ilask,2,-1 127 xa = yearlask(ilask) 128 xb = yearlask(ilask - 1) 129 if (xa <= Year .and. Year < xb) then 130 ! Obliquity 131 if (var_obl) then 132 ya = oblask(ilask) 133 yb = oblask(ilask - 1) 134 if (ya < yb) then ! Increasing -> max is the limitation 135 yc = max_obl 136 else ! Decreasing -> min is the limitation 137 yc = min_obl 138 endif 139 endif 140 141 ! Excentricity 142 if (var_ex) then 143 ya = exlask(ilask) 144 yb = exlask(ilask - 1) 145 if (ya < yb) then ! Increasing -> max is the limitation 146 yc = max_ex 147 else ! Decreasing -> min is the limitation 148 yc = min_ex 149 endif 150 endif 151 152 ! Lsp 153 if (var_lsp) then 154 ya = lsplask(ilask) 155 yb = lsplask(ilask - 1) 156 if (ya < yb) then ! Increasing -> max is the limitation 157 if (yb - ya > 300.) then ! If modulo is "crossed" then it is the other way around 158 !yb = yb - 360. 159 yc = min_lsp 160 else 161 yc = max_lsp 162 endif 163 else ! Decreasing -> min is the limitation 164 if (ya - yb > 300.) then ! If modulo is "crossed" then it is the other way around 165 !yb = yb + 360. 166 yc = max_lsp 167 else 168 yc = min_lsp 169 endif 170 endif 171 endif 172 iylask = ilask 173 exit ! The loop is left as soon as the right interval is found 174 endif 175 enddo 176 if (ilask == 1) then 177 write(*,*) 'The year does not match with Laskar data in the file ob_ex_lsp.asc.' 178 stop 149 179 endif 150 if(.not.var_ex) then 151 ex_not_found=.FALSE. 152 else 153 ex_not_found=.TRUE. 154 endif 155 if(.not.var_lsp) then 156 lsp_not_found=.FALSE. 157 else 158 lsp_not_found=.TRUE. 159 endif 160 161 max_obl_iter=999999 162 max_ex_iter =999999 163 max_lsp_iter=999999 164 165 !-------------------------------- 166 167 yc=max_obl 168 yc=min_obl 169 170 !-------------------------------- 171 172 do ilask=last_ilask,2,-1 173 xa=yearlask(ilask) 174 xb=yearlask(ilask-1) 175 ! Linear interpolation to find the maximum number of iteration 176 if((oblask(ilask-1).GT.max_obl).and. obl_not_found ) then 177 ya=oblask(ilask) 178 yb=oblask(ilask-1) 179 yc=max_obl 180 max_obl_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year 181 obl_not_found=.FALSE. 182 elseif((oblask(ilask-1).LT.min_obl).and. obl_not_found ) then 183 ya=oblask(ilask) 184 yb=oblask(ilask-1) 185 yc=min_obl 186 max_obl_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year 187 obl_not_found=.FALSE. 188 endif 189 if((exlask(ilask-1).GT.max_ex).and. ex_not_found ) then 190 ya=exlask(ilask) 191 yb=exlask(ilask-1) 192 yc=max_ex 193 max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year 194 ex_not_found=.FALSE. 195 elseif((exlask(ilask-1).LT.min_ex ).and. ex_not_found ) then 196 ya=exlask(ilask) 197 yb=exlask(ilask-1) 198 yc=min_ex 199 max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year 200 ex_not_found=.FALSE. 201 endif 202 if(max_lsp_modulo .and. lsplask(ilask).lt.180) lsplask(ilask)=lsplask(ilask)+360. 203 if(min_lsp_modulo .and. lsplask(ilask).gt.180) lsplask(ilask)=lsplask(ilask)-360. 204 if((lsplask(ilask-1).GT.max_lsp).and. lsp_not_found ) then 205 ya=lsplask(ilask) 206 yb=lsplask(ilask-1) 207 yc=max_lsp 208 max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year 209 lsp_not_found=.FALSE. 210 elseif((lsplask(ilask-1).LT.min_lsp ).and. lsp_not_found ) then 211 ya=lsplask(ilask) 212 yb=lsplask(ilask-1) 213 yc=min_lsp 214 max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year 215 lsp_not_found=.FALSE. 216 endif 217 enddo 218 219 print *, "Maximum obliquity accepted=", max_obl 220 print *, "Minimum obliquity accepted=", min_obl 221 print *, "Maximum number of iteration for the obl. parameter=", max_obl_iter 222 223 print *, "Maximum excentricity accepted=", max_ex 224 print *, "Minimum excentricity accepted=", min_ex 225 print *, "Maximum number of iteration for the ex. parameter=", max_ex_iter 226 227 print *, "Maximum lsp accepted=", max_lsp 228 print *, "Minimum lsp accepted=", min_lsp 229 print *, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter 230 231 year_iter_max=min(max_obl_iter,max_ex_iter,max_lsp_iter) 232 233 print *, "So the max. number of iteration (year) for the orbital parameter=", year_iter_max 180 181 ! Linear interpolation gives the maximum reachable year according to the limitation 182 ! Obliquity 183 if (var_obl) then 184 do ilask = iylask,2,-1 185 ya = oblask(ilask) 186 yb = oblask(ilask - 1) 187 if (ya <= yc .and. yc < yb) then 188 xa = yearlask(ilask) 189 xb = yearlask(ilask - 1) 190 max_obl_iter = int((max_obl - ya)*(xb - xa)/(yb - ya) + xa) - Year 191 exit ! The loop is left as soon as the right interval is found 192 endif 193 enddo 194 endif 195 ! Excentricity 196 if (var_ex) then 197 do ilask = iylask,2,-1 198 ya = exlask(ilask) 199 yb = exlask(ilask - 1) 200 if (ya <= yc .and. yc < yb) then 201 xa = yearlask(ilask) 202 xb = yearlask(ilask - 1) 203 max_ex_iter = int((max_ex - ya)*(xb - xa)/(yb - ya) + xa) - Year 204 exit ! The loop is left as soon as the right interval is found 205 endif 206 enddo 207 endif 208 ! Lsp 209 if (var_lsp) then 210 do ilask = iylask,2,-1 211 ya = lsplask(ilask) 212 yb = lsplask(ilask - 1) 213 if (ya <= yc .and. yc < yb) then 214 xa = yearlask(ilask) 215 xb = yearlask(ilask - 1) 216 max_lsp_iter = int((max_lsp - ya)*(xb - xa)/(yb - ya) + xa) - Year 217 exit ! The loop is left as soon as the right interval is found 218 endif 219 enddo 220 endif 221 222 print*, "Maximum number of iteration for the obl. parameter=", max_obl_iter 223 print*, "Maximum number of iteration for the ex. parameter=", max_ex_iter 224 print*, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter 225 226 year_iter_max = min(max_obl_iter,max_ex_iter,max_lsp_iter) 227 if (year_iter_max > 0) then 228 print*, "So the max. number of iteration (year) for the orbital parameter=", year_iter_max 229 else 230 write(*,*) 'The max. number of iteration (year) is not compatible with Laskar data in the file ob_ex_lsp.asc.' 231 stop 232 endif 234 233 235 234 END SUBROUTINE orbit_param_criterion -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3020 r3022 820 820 #endif 821 821 822 if (evol_orbit_pem) then822 if (evol_orbit_pem) then 823 823 call orbit_param_criterion(year_iter_max) 824 824 else 825 year_iter_max =Max_iter_pem825 year_iter_max = Max_iter_pem 826 826 endif 827 827 -
trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
r2985 r3022 4 4 ! Purpose: Recompute orbit parameters based on values by Laskar et al., 5 5 ! 6 ! Author: RV 6 ! Author: RV, JBC 7 7 !======================================================================= 8 8 IMPLICIT NONE … … 42 42 43 43 real :: Year ! Year of the simulation 44 real :: timeperi_ls! time of peri in ls44 real :: Lsp ! time of peri in ls 45 45 integer ilask ! Loop variable 46 46 real :: halfaxe ! Million of km 47 47 real :: unitastr 48 48 49 real xa,xb, xc,ya,yb49 real xa,xb,ya,yb 50 50 51 51 … … 61 61 #endif 62 62 63 Year=year_bp_ini+year_PEM+final_iter 64 65 timeperi_ls=timeperi*360/2/pi 63 Year = year_bp_ini + year_PEM + final_iter ! We compute the current year 64 Lsp = 360. - timeperi*360./(2.*pi) ! We convert in degree 66 65 67 66 ! Year=-953200 68 print 69 print 70 print 71 print 72 print *, "recomp_orb_param, Old lsp=", timeperi_ls73 print 67 print*, "recomp_orb_param, Old year=", year_bp_ini+year_PEM 68 print*, "recomp_orb_param, New year=", Year 69 print*, "recomp_orb_param, Old obl=", obliquit 70 print*, "recomp_orb_param, Old ex=", e_elips 71 print*, "recomp_orb_param, Old lsp=", Lsp 72 print*, "recomp_orb_param, Old timeperi=", timeperi 74 73 75 xc=Year 74 do ilask = last_ilask,2,-1 75 xa = yearlask(ilask) 76 xb = yearlask(ilask - 1) 77 if (xa <= Year .and. Year < xb) then 78 ! Obliquity 79 if (var_obl) then 80 ya = oblask(ilask) 81 yb = oblask(ilask - 1) 82 obliquit = (Year - xa)*(yb - ya)/(xb - xa) + ya 83 endif 76 84 77 do ilask=last_ilask,1,-1 85 ! Excentricity 86 if (var_ex) then 87 ya = exlask(ilask) 88 yb = exlask(ilask - 1) 89 e_elips = (Year - xa)*(yb - ya)/(xb - xa) + ya 90 endif 78 91 79 xa=yearlask(ilask) 80 xb=yearlask(ilask+1) 92 ! Lsp 93 if (var_lsp) then 94 ya = lsplask(ilask) 95 yb = lsplask(ilask - 1) 96 if (yb - ya > 300.) then ! If modulo is "crossed" 97 if (ya < yb) then ! Increasing 98 yb = yb - 360. 99 else ! Decreasing 100 yb = yb + 360. 101 endif 102 endif 103 Lsp = (Year - xa)*(yb - ya)/(xb - xa) + ya 104 Lsp = modulo(Lsp,360.) 105 endif 106 exit ! The loop is left as soon as the right interval is found 107 endif 108 enddo 109 halfaxe = 227.94 110 timeperi = Lsp*2.*pi/360. 111 periheli = halfaxe*(1 - e_elips) 112 aphelie = halfaxe*(1 + e_elips) 113 unitastr = 149.597927 114 p_elips = 0.5*(periheli + aphelie)*(1. - e_elips*e_elips)/unitastr 115 #ifndef CPP_STD 116 call call_dayperi(timeperi,e_elips,peri_day,year_day) 117 #endif 81 118 82 if(yearlask(ilask) .GT.Year) then 83 if(var_obl) then 84 ya=oblask(ilask) 85 yb=oblask(ilask+1) 86 obliquit=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb) 87 endif 88 if(var_ex) then 89 ya=exlask(ilask) 90 yb=exlask(ilask+1) 91 e_elips=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb) 92 endif 93 if(var_lsp) then 94 if(lsplask(ilask)-lsplask(ilask+1) .gt.200) then 95 if(lsplask(ilask).lt.lsplask(ilask+1)) then 96 lsplask(ilask)=lsplask(ilask)+360 97 else 98 lsplask(ilask+1)=lsplask(ilask+1)+360 99 endif 100 endif 101 102 ya=lsplask(ilask) 103 yb=lsplask(ilask+1) 104 timeperi_ls=ya*(xc-xb)/(xa-xb)+yb*(xa-xc)/(xa-xb) 105 timeperi_ls=modulo(timeperi_ls,360.) 106 halfaxe=227.94 107 timeperi=timeperi_ls*2*pi/360 108 periheli = halfaxe*(1-e_elips) 109 aphelie = halfaxe*(1+e_elips) 110 unitastr=149.597927 111 p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr 112 #ifndef CPP_STD 113 call call_dayperi(timeperi,e_elips,peri_day,year_day) 114 #endif 115 endif 116 exit 117 endif 118 enddo 119 120 print *, "recomp_orb_param, Final year of the PEM run:", Year 121 print *, "recomp_orb_param, New obl=", obliquit 122 print *, "recomp_orb_param, New ex=", e_elips 123 print *, "recomp_orb_param, New timeperi=", timeperi 119 print*, "recomp_orb_param, Final year of the PEM run:", Year 120 print*, "recomp_orb_param, New obl=", obliquit 121 print*, "recomp_orb_param, New ex=", e_elips 122 print*, "recomp_orb_param, New timeperi=", timeperi 124 123 125 124 END SUBROUTINE recomp_orb_param
Note: See TracChangeset
for help on using the changeset viewer.
