MODULE orbit_param_criterion_mod !======================================================================= ! ! Purpose: Compute the maximum nimber of iteration the PEM can do based ! on the stopping criterion given by the orbital parameters ! ! Author: RV !======================================================================= IMPLICIT NONE CONTAINS SUBROUTINE orbit_param_criterion(year_iter_max) #ifdef CPP_IOIPSL use IOIPSL, only: getin #else ! if not using IOIPSL, we still need to use (a local version of) getin use ioipsl_getincom, only: getin #endif 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) ! Number of line in Laskar file 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) logical max_lsp_modulo,min_lsp_modulo ! Variable to check if we are close to the 360 degree modulo for lsp parameter real xa,xb,ya,yb,yc ! ********************************************************************** ! 0. Initializations ! ********************************************************************** Year=year_bp_ini+year_PEM ! We compute the current year timeperi_ls=timeperi*360/2/pi ! We convert in degree call ini_lask_param_mod(nlask) ! Allocation of variables 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 ! We read the file 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.5 max_change_ex=0.1 max_change_lsp=40. CALL getin('max_change_obl', max_change_obl) CALL getin('max_change_ex', max_change_ex) CALL getin('max_change_lsp', max_change_lsp) 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 max_lsp_modulo=.false. min_lsp_modulo=.false. if(max_lsp.gt.360) max_lsp_modulo=.true. if(min_lsp.lt.0) min_lsp_modulo=.true. !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=999999 max_ex_iter =999999 max_lsp_iter=999999 !-------------------------------- yc=max_obl yc=min_obl !-------------------------------- do ilask=last_ilask,1,-1 xa=yearlask(ilask) xb=yearlask(ilask-1) ! Linear interpolation to find the maximum number of iteration if((oblask(ilask-1).GT.max_obl).and. obl_not_found ) then ya=oblask(ilask) yb=oblask(ilask-1) yc=max_obl max_obl_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year obl_not_found=.FALSE. elseif((oblask(ilask-1).LT.min_obl).and. obl_not_found ) then ya=oblask(ilask) yb=oblask(ilask-1) yc=min_obl max_obl_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year obl_not_found=.FALSE. endif if((exlask(ilask-1).GT.max_ex).and. ex_not_found ) then ya=exlask(ilask) yb=exlask(ilask-1) yc=max_ex max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year ex_not_found=.FALSE. elseif((exlask(ilask-1).LT.min_ex ).and. ex_not_found ) then ya=exlask(ilask) yb=exlask(ilask-1) yc=min_ex max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year ex_not_found=.FALSE. endif if(max_lsp_modulo .and. lsplask(ilask).lt.180) lsplask(ilask)=lsplask(ilask)+360. if(min_lsp_modulo .and. lsplask(ilask).gt.180) lsplask(ilask)=lsplask(ilask)-360. if((lsplask(ilask-1).GT.max_lsp).and. lsp_not_found ) then ya=lsplask(ilask) yb=lsplask(ilask-1) yc=max_lsp max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year lsp_not_found=.FALSE. elseif((lsplask(ilask-1).LT.min_lsp ).and. lsp_not_found ) then ya=lsplask(ilask) yb=lsplask(ilask-1) yc=min_lsp max_ex_iter=xa*(yc-yb)/(ya-yb)+xb*(ya-yc)/(ya-yb)-Year 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_lsp print *, "Minimum lsp accepted=", min_lsp 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