MODULE orbit_param_criterion_mod IMPLICIT NONE CONTAINS SUBROUTINE orbit_param_criterion(year_iter_max) USE temps_mod_evol, ONLY: year_bp_ini, year_PEM, var_obl, var_ex, var_lsp #ifndef CPP_STD USE planete_h, ONLY: e_elips, obliquit, timeperi #else use planete_mod, only: e_elips, obliquit, timeperi #endif USE comconst_mod, ONLY: pi USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, & ini_lask_param_mod,last_ilask IMPLICIT NONE !-------------------------------------------------------- ! Input Variables !-------------------------------------------------------- !-------------------------------------------------------- ! Output Variables !-------------------------------------------------------- integer,intent(out) :: year_iter_max ! Maximum number of iteration of the PEM before orb changes too much !-------------------------------------------------------- ! Local variables !-------------------------------------------------------- real :: Year ! Year of the simulation real :: timeperi_ls ! Year of the simulation integer nlask,ilask!,last_ilask !Loop variable parameter (nlask = 20001) real max_change_obl,max_change_ex,max_change_lsp ! Percentage of change that is considered to be acceptible real max_obl,max_ex,max_lsp !Maximum value of orbit param given the acceptable percentage real min_obl,min_ex,min_lsp !Maximum value of orbit param given the acceptable percentage real max_obl_iter,max_ex_iter,max_lsp_iter !Maximum year iteration before reaching an unacceptable value logical obl_not_found, ex_not_found,lsp_not_found !Loop variable (first call) ! ********************************************************************** ! 0. Initializations ! ********************************************************************** Year=year_bp_ini+year_PEM timeperi_ls=timeperi*360/2/pi call ini_lask_param_mod(nlask) print *, "orbit_param_criterion, Year in pem.def=", year_bp_ini print *, "orbit_param_criterion, Year in the startpem.nc =", year_PEM print *, "orbit_param_criterion, Current year=", Year print *, "orbit_param_criterion, Current obl=", obliquit print *, "orbit_param_criterion, Current ex=", e_elips print *, "orbit_param_criterion, Current lsp=", timeperi_ls open(73,file='ob_ex_lsp.asc') do ilask=1,nlask read(73,*) yearlask(ilask),oblask(ilask), & exlask(ilask),lsplask(ilask) yearlask(ilask)=yearlask(ilask)*1000 if(yearlask(ilask).GT.Year) then last_ilask=ilask+1 endif end do close(73) print *, "Coresponding line in the ob_ex_lsp.asc file=", last_ilask ! 5% max change case max_change_obl=0.05 max_change_ex=0.05 max_change_lsp=0.05 max_obl=obliquit*(1+max_change_obl) min_obl=obliquit*(1-max_change_obl) max_ex=e_elips*(1+max_change_ex) min_ex=e_elips*(1-max_change_ex) max_lsp=timeperi_ls*(1+max_change_lsp) min_lsp=timeperi_ls*(1-max_change_lsp) !End of 5% max change case !Constant max change case max_change_obl=0.1 max_change_ex=0.1 max_change_lsp=40. max_obl=obliquit+max_change_obl min_obl=obliquit-max_change_obl max_ex=e_elips+max_change_ex min_ex=e_elips-max_change_ex max_lsp=timeperi_ls+max_change_lsp min_lsp=timeperi_ls-max_change_lsp !End Constant max change case ! If we do not want some orb parameter to change, they should not be a stopping criterion, ! So the number of iteration corresponding is set to maximum if(.not.var_obl) then obl_not_found=.FALSE. else obl_not_found=.TRUE. endif if(.not.var_ex) then ex_not_found=.FALSE. else ex_not_found=.TRUE. endif if(.not.var_lsp) then lsp_not_found=.FALSE. else lsp_not_found=.TRUE. endif max_obl_iter=999999999999 max_ex_iter =999999999999 max_lsp_iter=999999999999 do ilask=last_ilask+1,1,-1 if((oblask(ilask).GT.max_obl).and. obl_not_found ) then max_obl_iter=(max_obl-oblask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & / (oblask(ilask-1)-oblask(ilask)) obl_not_found=.FALSE. elseif((oblask(ilask).LT.min_obl).and. obl_not_found ) then max_obl_iter=(min_obl-oblask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & / (oblask(ilask-1)-oblask(ilask)) obl_not_found=.FALSE. endif if((exlask(ilask).GT.max_ex).and. ex_not_found ) then max_ex_iter=(max_ex-exlask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & / (exlask(ilask-1)-exlask(ilask)) ex_not_found=.FALSE. elseif((exlask(ilask).LT.min_ex ).and. ex_not_found ) then max_ex_iter=(min_ex-exlask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & / (exlask(ilask-1)-exlask(ilask)) ex_not_found=.FALSE. endif if((lsplask(ilask).GT.max_lsp).and. lsp_not_found ) then max_lsp_iter=(max_lsp-lsplask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & / (lsplask(ilask-1)-lsplask(ilask)) lsp_not_found=.FALSE. elseif((lsplask(ilask).LT.min_lsp ).and. lsp_not_found ) then max_lsp_iter=(min_lsp-lsplask(ilask)) * (yearlask(ilask-1)-yearlask(ilask)) & / (lsplask(ilask-1)-lsplask(ilask)) lsp_not_found=.FALSE. endif enddo print *, "Maximum obliquity accepted=", max_obl print *, "Minimum obliquity accepted=", min_obl print *, "Maximum number of iteration for the obl. parameter=", max_obl_iter print *, "Maximum excentricity accepted=", max_ex print *, "Minimum excentricity accepted=", min_ex print *, "Maximum number of iteration for the ex. parameter=", max_ex_iter print *, "Maximum lsp accepted=", max_obl print *, "Minimum lsp accepted=", min_obl print *, "Maximum number of iteration for the lsp. parameter=", max_lsp_iter year_iter_max=min(max_obl_iter,max_ex_iter,max_lsp_iter) print *, "So the max. number of iteration (year) for the orbital parameter=", year_iter_max END SUBROUTINE orbit_param_criterion !******************************************************************************** END MODULE orbit_param_criterion_mod